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ABS TRACT 


This thesis treats the subject of digital filters in two 
parts. First the analytical tools and synthesis techniques 
which are used for digital signal processors are explained 
and illustrated with examples. Emphasis is placed on 1) 
the generation and use of the digital transfer function, 2) 
design procedures in the frequency domain using methods of 
continuous-filter design as intermediate steps, and 3) 
real-time digital-filter implementation schemes. 

In the second part of the thesis a machine-language 
computer program 1S presented for the real-time simulation 
of digital filters on a hybrid computer. The program is 
described in detail and experimental results for several 
filter designs and realization schemes are reported. It is 
believed that such a computer program, with its inherent 
ability to accurately simulate a special-purpose computer and 
to provide external control of word length and clock frequen- 
cy, could be used to experimentally determine specifications 


for the LSI implementation of digital.filters. 
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Ll. -ENGRODUCTION 


The process of signal waveform filtering in linear sys- 
tems may generally be defined as the application of a linear 
operator to an input waveform as it passes through a system. 


The transformation may occur in the time domain such as 


Si 4 


eels. oy extrapolation, differentiation and integra- 
Pon: or it may occur in the frequency domain such as low- 
pass filtering or smoothing, bandpass filtering, spectrum 
shaping and spectral decomposition. Operations such as 

these on continuous-time signals using analog filters are 
defined ee by linear constant-coefficient differ- 
ential equations. Analysis of these operations is made 
possible in the s domain by the LaPlace transform, and in the 
12 domain by- the Fourier transform. Similarly linear trans- 
formation#®on discrete-time Signals may be defined by linear 


constant-coefficient difference equations, and discrete-time 


filters may:-be analyzed and synthesized using the z-transform. 


AeA GENERAL DISCUSSION OF DIGITAL FILTERS 

A digital: filter 1s defined as a time-invariant linear 
operator on discrete-time signals. This linear operator is 
basically a computational process by which a sampled signal 
or sequence of numbers, called the input, is transformed into 
a second sequence of numbers, called the output. Occasionally 
the term "digital filter" includes the means for performing 


the analog-to-digital and digital-to-analog conversions. 


te: 


dgll 


However, since some signals are already digital in nature, 
"digital filter" ines 2 11 refer only to the compu- 
tational process described above. 

Digital filters are time-invariant in the sense that they 
are described by linear constant-coefficirent ditrerence 
equations. Any given set of coefficients uniquely specifies 
a filter, Je any change in these coefficients would neces- 
Sarily result in a change in the specified filter. MThus, a 
digital filter is time-invariant. 

TE cen eT On of the computational processes required in 
digital filters is accomplished using digital components such 
as shift registers, adders and gates. Several fundamental 
advantages to be gained using digital hardware are reduction 
in cost, higher speed, and increased reliability. Digital 
Signal-processing techniques offer the following additional 
advantages over analog signal processing: 


1. Very predictable, drift-free performance of arbitrarily 
high precision can be obtained. This permits the 
realization of stable filters with very high Q's. 


Dus Theré are no impedance-matching problems since digital 
circuits are inherently isolated from each other. 


3. There is no restriction on critical filter frequencies. 
In particular, filters for frequencies of fractions of 
hertz are easy to construct. 


4. A greater variety of filter functions is possible since 
implementation is essentially a matter of numerical 
coefficients vice inductors and capacitors. 


5. Filter response can be changed easily by varying the 
proper coefficients. 


6. There is an intrinsic possibility of time-sharing the 
major filter sections. For example, a many-Section 
filter can be implemented by building only one section 
and multiplexing that section by time-switching it. 

The filter coefficients can be changed for each section. 


F Alternatively, a given filter may operate on more than 
one input sequence by time multiplexing. 


These advantages in conjunction with the potential impact 
of large-scale circuit integration (LSI) on real-time digital 
filters assures wide application for these devices in such 
areas as audio systems, speech processing, seismology, ocean- 


ography, radar systems and data communications. 


B. HISTORICAL BACKGROUND 

The processing of discrete signals using linear filters 
or weighting sequences is traced back to the 1600's by 
Kaiser [1]. He discusses the early use of classical rasa ONL 
methods amd the subsequent evolution of the z-transform cal- 
culus. The theory of the z-transform and its application to 
the early digital signal processors are contained in the 
extensive literature on sampled-data control systems published 
in the 1950's. See, for example, Refs. 1, 2, and 3. 

More recently the general-purpose digital computer became 
a vehicle for the simulation of continuous filters, and the 
development of digital-filtering techniques per se was greatly 
accelerated.- As computer technology advanced, particularly in 
the areas of memory size and computational speed, the real- 


Ayate 


time implementation of digital filters on general-purpose 


od 


‘computers became possible.? Now, with the advent of 


vat 


asa 


‘ Indeed, the full impact of high-speed convolution via 
the fast Fourier transform (FFT) on real-time digital fil- 
tering 1s not yet known. | 


~- 
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large-scale circuit integration, the realization of digital 
filters in special-purpose hardware is becoming practices 
and economical. 

In short, the prospect of building real-time digital 
filters via LSI has generated considerable interest in the 
general field of digital Signal processing. References 4 - 8 
contain much of the most recently published material on this 


relatively new subject. 
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II. MATHEMATICAL CONSIDERATIONS 


The analysis and design of digital filters requires a 
moderate degree of familiarization with the mathematics of 
discrete-time signals. In this section the the@ry of the 
sampling process and the z-transform calculus is reviewed, 
and their relation to the digital processing of signals is 
presented. Particular attention is paid to the methods of 
ee ewating the digital transfer function, H(z); because, as 
will be shown later, this operation comprises the essence of 


the digital-filter design problem. 


Pe THE SAMPLING PROCESS 

The operation of sampling a continuous-time Signal can 
be thought of as a form of impulse modulation. If a signal 
£(t) is sampled every T seconds, then the sampling rate is 
W. = 27/T, and the sampled signal is a train of equally 


Ss 


spaced impulses with varying amplitude given by 


£*(t) = £ £(nT)6(t-nt). (2a) 
n= Gene 


In the analog-to-digital converter the amplitude of f(t) 
at each sample time, £(nT), 1s “measured," and a digital 
word is then formed by encoding a sequence of binary digits 
corresponding to that eo value. It is immediately seen 
that the accutacy of the "measurement" will be limited by the 

length of the computer word which the digital filter is 
Meeucned to handle. This phenomenon is known as quantization 


(4, 5], and its effects are discussed in Section IV, D. 


ys 


Returning to Eq. (2-1), it 1S evadenitemghat, in the time 
domain, the sampled function £*(t) equals the original 
Pinciteonemulcaplied by sche modulating muneELOn, o(t=nl)s 
Consequently, in the frequency domain, the spectrum of f£¥*(t) 
will be equal to the convolution of an infinite number of 
equally spaced spectral lines with the original signal spec- 
trum. Hence, the periodic spectrum shown in Fig. 2-1 is 
obtained. It is important to note that in a real sampler the 
modulating pulses are not mathematical delta functions, but 
pulses with a finite width. This explains the decreasing 
amplitude in the spectra along the frequency axis on either 
side of zero frequency. The Shape of this attenuating 
function is the familiar sin x/x shape, and its derivation 
can be formalized using the Fourier transform. 

Of particular significance is the phenomenon of fre- 
quency folding or aliasing. Again refer to Fig. 2-l. If 


Ww 1s the highest frequency component of the original signal 


ba eel 


spectrum and w. < 2W ar distortion will appear in the sampled 


Signal spectrum because of the overlapping of signal compo- 


nents. In general, therefore, 1n order to recover the 


original spectrum with a low-pass filter, that spectrum must 
| ae 


See 





- 


be bandlimited to twi/2. This is known as Nyquist's Sampling 
Theorem [1-3] and 20, is called the Nyquist frequency (minimum 
sampling aceon In practice there will always be some 
degree of frequency folding. Hence the total amount of in- 


FOLrmMatien ts never recoverable. 
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It should be emphasized at this, point that the sampling 
frequency must be the clock fin ca ey Lerethesdigita | £1 Item. 
Pins teases aaa sel) numerical computation that is required to gen- 
erate one output value must be completed in one sampling 
period. Hencer the computation time of the digital £1125 


sets an upper limit on the sampling frequency. 


B. LINEAR DIFFERENCE EQUATIONS 

If an input sequence is denoted by a and ag output 
sequence by va then the digital filtering of x to produce 
Y, May be expressed by linear constant-coefficient difference 


equations of the form 
eee + = eee 
Yn * D5 ae a PuYn—-M sa Ty eee, u T ay*n-n’ (2-2) 


Each term in the above equation is a sequence of numbers. 
The sequence iar differs’ fromey_ only in=that it is delayed 
enevunet (T)yam times Bq@uation (2-2), written as a funetion 


of T, becomes 
| (2538 ) 
y Cri) ae biy(nt-T) . cee + yy (nT-MT) = a,x (nT) a,x(nT-T) 


+ eoe 4 ax (nT-NT) . 
Equations (2-2) and (2-3) are equivalent. The notation of 
Eq. (2-2) is often preferred because of simplicity. The 
interpretation of the generalized difference equation is as 
follows:. at time t=nT the new value in the output sequence 
Yn, can be determined from the current input sample, the most 
recent value in the input sequence ar and a linear combina- 


an 


tion of the past input and output values which are available 
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for computation. The constant coefficients, a, and b., are 


determined by the particular filtering function desired. The 


ea cen TEED 


specification of these coefficients becomes the essence of 
the digital-filter design problem. 


An example of a first-order linear difference equation is 
y(nT) = b,y(nT-T) + X(t); (2-4) 


(Here first-order refers to the fact that the highest order 


delay in the equation is one.) If x(t) is the complex 
| He Accept, for the moment, 
that y(nT-T) can be expressed as y(nt)e 2°, ves, ymT-T) 


sinusoid e Give Ht Aree SS 
is y(nT) delayed one sampling period in time. (This idea 
will be formalized in the next section.) Then Eq. (2-4) can 


be written as 


jnwT JOT 
7 e _ x(nT)e he 
1 - b,e e ~ - 


Equation (2-5) shows that, in the steady state, the output 
is also a complex sinusoid, but modified by a "transfer 


Sonat) . 


fM@nction." The transfer function can be defined as H(le 


Its magnitude is 


13 ik 
pH (eI?) | : ; os (2-6) 
(1 + by - b,coswT) 
ud 1k 
and the phase is 
. 2 mall SinwT = 
¥ = er ~ tan“ (Soswr=b:) ° (2p) 


{ 


Equation (2-6) illustrates the frequency-selective properties 


of the transfer function, Hie. Note that the magnitude 
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response 1S periodic in frequency with a period equal to 
0) l= 27™/T. These results are, in fact, quite general and 
the concept of a discrete-signal transfer function will be 


discussed in detail in the next section. 


C. THE Z-TRANSFORM 

It bas been stated that digital filtering involves the 
computational process of obtaining a solution to linear 
constant-coefficient difference equations of the form given 


in Eq. (2.2). Furthermore, it was implied that the steady- 


— 





state solution to Eq. (2-2) for a Sinusoidal input can be 
obtained by applying a transfer function to the input signal. 
Transformations which permit the formulation of .a discrete- 
time signal transfer function are called sampled-data or 
ZeemeansrOrms. “The temm “digmtal transter function” wile 
used to describe a discrete-time transfer function of the 
complex variable z. 


1. Definition of the Z-Transform 


Given a discrete signal!f, Yepresented by a Sequamce 


of numbers f,; its z-transform, F(z), is defined by the power 
series 
a r § ar 
Fae = » {2 ’ (2735 
n=o 
where z is interpreted as a complex transform variable. The 


Variable z plays a role similar to that of the LaPlace 
transform variable s. The operation of taking the z-transform 
of a sequence is denoted by 


Paez) = Z(f_). (2-9) 
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As a simple example, if t. were the constant signal 
fe = Hie n = O, Alb 25 ok r (2-10) 


then the z-transform is found from Eq. (2-8) to be the 
geometric series 


Z(1l) = ~~ = | Zale || (2-11) 


n=o | 2 ee 
The theory of the z-transform 1s well documented 
{2,3]; detailed proofs and further examples will not be given 
here. Rathew, an important theorem and several conclusions 
which apply 46 digital filters are the following: 


1. The z-transform of a delayed sequence can be derived 
Leone derantcicw. Phe result is 


ae) See ene ei 


—_ 


2. The z-transform of the real part of a complex sinusoid 
is 
sa ces b za 
Z(acos nb) = SS ’ |z|>a ° (2-13) 
(l-ae?"z ~) (l-ae ?~z ~) 


3. The z-transform of the imaginary part of a complex 
Sinusoid is 


| =i 
g(a"sin nb) = ——2 Sin bz 9-414) 


(1-aediae. |) (12a ae 


Note that the z-transforms of sinusoids yield poles 
in complex-conjugate pairs at a radius depending on 
the damping factor. j 


4. The z-transform of any function of the general form 


f = Olen Sin(nb + c) (2—-Fay 


and linear combinations_of these will always be a 
rational function of z™~~, and such functions may be 
broken down into terms of the form of Eqs. (2-13) "ama 


(2-14). This class of signals corresponds to the class 
of continuous-time signals with rational LaPlace 
transforms. 


Pagll 


The correspondence between the z-transform and the 
LaPlace transform will now be shown. From Eq. (2-1) a 
sampled signal is 


oe i 2 f(nT)d(t=ng)- (2—Tey 
fio 


Taking the Laplace transform of both sides of Eq. (2-16) 


yields 


Tx = z = 
F¥*¥(s) nes (2-17) 


2 Af 
This result rts identical in form with the defini trom o mene 


z-transform, Eq. (2-8), provided the substitution, 


a S72 Riis) 


is made. The complex variable z+ May be interpreted as 
a "unit delay operator" in Eq. (2-8); or, ina different 
sense, aS an ordering variable whose exponent represents the 
position of a pulse in a sequence. This latter interpreta- 
tion, sometimes referred to as a’ "generating function" [2], 
will become Bhat after the discussion of the z-transform 
inverse. 

2>— ThewmDigiital Trangicnaluneduens 


The output of a digital filter is again written as 


ee OP gee (ee Oe a elo ene 


? ee AN nen 


(Zao) 


Taking the z-transform of each term and applying the 


delayed sequence theorem gives 


ONT! 


M 


= - 2 = 
Z(y_) + b,z Z(y.) + ess + DiyZ Z(y,) = ao2 (x) + az Z (x) 
= 
+r ate Anz Z(x,). 
Factoring the common terms and simplifying Eq. (2-20), the 
result 1s N . 
y aa 
i=o 7 
Y(2) _ psec. OO _ 
ma - H(z)- (221) 
dg IS lee! oy 72 
ii 
as 


H(z) is defined as the ratio of the z-transformed output to 
the z-transformed input and will generally be called a 
digital transfer function. 

Aff ther approach to the formulation of a digital 


transfer function is one using the notation of a linear 


operator H. If an output sequence is given by 


Yu = 48%, + 8y¥%n-7 7 °° CZ=2 2) 


which is simply the input sequence multiplied by a weighting 


function, then the notation can be simplified to 


Ye ‘= H(x,) S18 n-1, iil a (2-23) 


where H defines the filter. If h, represents the impulse 
response of one filter, (1g erage oe response to an input 
sequence which is unity at n=0 and is ZeLO at all other 
times), then the time response ee obtained from the 


discrete convolution theorem [4], 
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y(nT) = ys 9x (kaye (nT=-kT) (2-24) 
=O 


or, equivalently, 
y = a x AL iy A (Z—Z) 


Taking the z-transform of both sides of Eq. (2-25) and again 


using the delayed-sequence theorem, 


ity.) == x4 (h,) ts xz °° 2(h_) +o see | 
or 


eZ) = H(z) a 
k=o 


bec CeO eee (2-26) 


Again, H(z) is defined as the digital transfer function of 
the dagieta)] *faiiter.H., 

3. The Z-Transform Inverse 

The transform of the output sequence Y(z) is seen to 

be obtained from the relation Y(z)=H(z)X(z). Assuming that 
both H(z) and X(z) are known and are rational polynomials in 
z, there remains the problem of determining the sequence yy 
from Y(z). Generally there are three methods for determining 
the inverse of a z-transform [2,3]. The first two are the 
use Of a z-transform inversion formula and the partial frac- 
tien expansion. The inversion formula is the most general, 
the partial fraction expansion requiring factorization of 
the denominator polynomial into terms for which transform 
pairs are known. 

A third method, and usually the easiest, is the power 
series or long division method. Consider the function Y(z) 


-1 
written in terms of Zz as 
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YZ) = 2 __+_—_+_.—__—- : (Ze) 
Do te bz ch b3z pues 


Tf the numerator is divided by the denominator using long 


division, the result is 
Y (2) ee eraser q,2 + q2z + eee , (2-28) 


This power series gives the output directly. 


pate) =e COME q, (t-T) + qo (t-2T) fueueees) =e (2-29) 


Here the variable z is interpreted as the output Sequence 
"generating function" mentioned earlier. Such a routine is 


easily performed on a computer or calculating machine. 


4. Sinusoidal-Steady-State Solution of Lanear Difference 
Equations 


In continuous-time linear systems the labor of taking 
inverse LaPlace or Fourier transforms can be avoided for the 
Special case of sinusoidal input functions. The output of 
Such a system in the steady state will always be a sinusoid 
of the same frequency, but whose magnitude and phase are 
modified by the transfer function of the system. The solu- 
tion for the magnitude of the Output over a specified fre- 
quency range is called the frequency response of the system. 
This function is found by evaluating the system transfer 
function H(s) or F(jw) along the real frequency axis. It 
will now be shown that the response of a discrete-time linear 
system can be found by evaluating the digital transfer func- 


tion H(z) around the unit circle in the z-plane. Hence, 


a) 


for the very common case of an input sequence being a 
uniformly sampled sine wave, the output sequence may be 
determined without calculating a z-transform inverse. 

In Section Il, B Eq. (2-5) showe@® Gear the eutpucter 
a discrete first-order system. to a sinusoidal input was the 
same Sinusoid modified by a frequency-sensitive transfer 
jwoT 


). Applying the substitution z=eot, with 


function H(e 
S=j]W, permits the transfer function to be written as a 


Pine ONO fr a2. 


gE 
neue) _ ae = ee = Tl ez). (2-30) 
ais i 
e - bo 1 - bz 


If the original difference equation, Eq. (2-4), had been 
solved using the z-transform, the resultant digital transfer 
function would have been identical to the one in Eq. (2-30) 
above. These results are not just coincidental, but are 
required by the definition of the z-transform. Therefore, it 
can be stated generally that the frequency response of a 
onKe pier. me filter can be and by evaluating the digital trans- 
fer function H(z) of that filter at telcos 

This discussion is illustrated graphically in Fig. 2-2. 
The mapping of the s-plane into the z-plane 1s accomplished 
by substituting S= OF Te Oss ECiele2 ald eee eT 


(Og a) = oot jot 


Zo =e -e ‘ (2-31) 


2In this context the term digital filter is intended to 
refer only to frequency-selective filters in the analog 
sense, vice the general time-domain filters mentioned in 
Sectrem Is AX 
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Figure 2-2. S-Plane to Z-Plane Transformation 


Za 


where T=27/0 For wi sw /2 the right portion of the s- 
plane (positive 0) maps into the area outside the unit circle 
in the z-plane. Similarly, the left portion of the s-plane 
(negative ¢) maps into the area inside the unit circle. Just 
as the frequency response of a continuous-time filter is 
determined by the value of its transfer function on the jw 
axis, the frequency response of a discrete-time filter is 
determined by the value of its transfer function on the unit 
circle. If the poles and zeros of a digital transfer func- 
tion are plotted in the z-plane,> then the magnitude of the 
response at the frequency @, in Fig. 2-2 18S li/L2L3 andere 
phase 1S G1-02-0;. These results are obtained using conven- 
tional pole-zero graphical techniques. It should be noted 
that the magnitude response repeats after w reaches Wor 2Wo, 
etc., and 1S symmetrical around NW. This result agrees with 
the periodic spectrum of the sampled signal discussed previ- 
ously. In general the phase characteristic may or may not 
repeat {9}. Also, the stability of a digital filter requires 
that the poles of the digital transfer function be located 
within the unit circle in the z-plane. 

The frequency response of a digital transfer function can 


I 


eae ar ee 


be evaluated analytically by substituting e 
the functional expression H(z). Assuming H(z) can be factored 


into products of second-order expressions, then 


Poles angurecos Of a duairat etranster function are 
found, as expected, by the factorization of the rational 
polynomials in Zz. 


A +A a e )-a 
O 1 


1 + Bie 


H(e7J47) Y = 
JWT ie EyOu an 


] WT =r 
a Ad + Ae 


= 5 je 1 
aF By =p Boe 


Expanding the exponential terms and regrouping yields 


Ae 
Oo 


eJut 


[(A,+A,) coswT+A,]+j[(A,-A Soa 


9) 


@7J4T) 
[ (B,+1) coswI+B,) ]+j[(1-B,) sinwT] 


H ( ° (2-32) 


Computer Program 1 may be used to evaluate the magnitude of 
the right side of Eq. (2-32). The program is written in the 


FORTRAN IV language. 
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D. S-PLANE TO Z-PLANE TRANS FORMATIONS \” 

Many of the digital-filter design techniques to be dis- 
cussed in Section III involve the digitalization of analog 
filters. These indirect design methods require transforma- 
tions from functions of s to functions of z. This section 
discusses the mechanics of various mapping transformations 
from the s-plane to the z-plane and compares the relative 
merits of each transformation. The formulas given for compu- 
tation of coefficients are not derived, but reproduced from 
the literature [6]. 

i. ‘The otandard Zeneca 

This transformation is simply the extension of the 
definition of the z-transform. A proper rational function 


of s 1S expanded into partial fractions and the z-transform 


of each term is calculated. If RR represents the real part 
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of a residue in the partial fraction expansion, and RI 
represents the imaginary part, then real terms are trans-— 


formed according to 





A 
= € ——~7 ’ (25537 
S ae 1+ Byz 
it 
where 
A = T (RR) 
O 
By =. = exp lull). 


Complex conjugate terms may be combined and transformed 
according to the transform pair 


A, + Aa 
Sey (2-34) 


-l -2 
1 + Biz + Boz 


RR + jRI _ RR - jRI 
Su - JV )ClCUS SU 


where 


= 2T(RR) 
= 2T exp(uT) [RR cos(vT)+RI sin(vT) ] 
2 Queso! UT igechda7T } 


OD WwW Pp Pp 
SE ee O 


= exp(2uT). 

Since the transformation involves the relation, zt =e ; 
the mapping of the s-plane into the z-plane takes place 
exactly as was previously described for Fig. 2-2. Notice 
that the transformation maps each strip in the s-plane, 
bounded by w = (n+1/2) 0, and w = (n-l/2)w., where n is an 
integer, into the entire z-plane. Thus, the mapping of each 
Succeeding strip is superimposed on the baseband mapping func- 
tion. If the transfer function H(s) has any significant 
response above w = Wf 2 (1.e., 1f£ the product of the zero 
Wecters Givided by the product of the pole vectors is not 


negligible for ww /2, frequency folding will occur and the 


2118) 


transformation may be unsatisfactory. Consequently the 
standard z-transform is normally used only for bandlimited 
functions. The problem is related to the one discussed for 
the sampled signal. In the present case, however, the 
transfer function itself must be bandlimited in addition to 
the input signal. In cases where H(z) is not sufficiently 
bandlimited a low-pass "guard" filter may be inserted in 
cascade to permit a successful transformation. 
2. The Bilinear Z-Transform 

To circumvent the problem of frequency folding a 
transformation is sought which will map the entire s-plane 
into the z-plane. As a first step, the entire s-plane is 
mapped into an s'-plane bounded by the lines s'‘' = +Jwi/2. 
The mapping relation is 

s = = tanh (£4). (2-35) 

Note that this transformation also maps the entire s-plane 
into each succeeding strip bounded by w' = at 3) and 
wi = Mss) w” for integer values of n. This condition is 
necessary Since the subsequent function of z must be periodic 
in frequency. Then, substituting z = aes Inver qs (24655 , 


the bilinear z-transform becomes 


ay Sa J -l 
Ss = = tanh (8.7 2 tine __) = tine ; (2286 
2 T 1+te ) 1 (lars a 


This algebraic transformation does map the entire left half 
of the s-plane into the interior of the unit circle in the 
z-~plane. Folding errors are eliminated since no folding 


occurs. The price paid for this feature is the nonlinear 


yall 


warping of the frequency scale. If s' = jw’ and s = jw in 
Eqs@ @62=35)eecthem 


TY gable 
a) = Can 7) ° (Pm oF) 


Since w' is the frequency which is really mapped into the 
Z-plane, the warping of the frequency scale according to the 
tangent function results. However, satisfactory results can 
usually be obtained by prewarping the design frequency before 
performing the transformation. 

3. The Matched Z-Transform 

This transform generates a digital transfer function 

with poles and zeros matched to those of the continuous 


function. The mapping transformation 1s 


o_ > i a (2-38) 
Then real poles or zeros transform according to 


-l 
ae 1 + A,z ‘ (2-39) 


where 
A, = - exp(uT) 
Complex poles or zeros are combined into second-order factors 


and transformed according to 


(s-u)* # v2 + Lt Ae + aie 7, (2-40) 
where 
Ay = -2exp(uT) cos vT) 
A; = exp (2uT) 
The matched z-transform preserves the shape of the 
frequence response characteristics. Golden [6] points out 


that it may sometimes be necessary to multiply the resultant 
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Mey by terms of the form (1 4 27-)" 


to achieve satisfactory 
results. This operation is equivalent to inserting n addi- 
tional zeros at the half-sampling frequency into the digital 
transfer function. 
4, Other Transformations 
Kaiser [1] discusses various other s-plane to z-plane 
transformations which may be used to achieve specific results. 


One particular transformation of interest is 


2 - 22 cosw oT + J] 


so ~~ (2-41) 
"Hii 


~~ 


This transformation, which is a variation of the bilinear 
z-transform, maps the entire imaginary axis of the s-plane 
‘onto both the top and bottom arcs of the unit circle in the 
z-plane. Setting s = 0 and solving the numerator of (2-41) 
for z shows that the origin of the s-plane is mapped into 
the points e + jot. Thus, applying this transformation to 
a low-pass H(s) yields a digital bandpass H(z) with center 


frequency of 0 If we denote the analog cutoff frequency 


by sn and the digital cutoff frequency by Wry then (2-41) 
can be written in equation form as 
20° 7 JW p,Tt 
% oO” -2e i cosw T + 1 
JWn aa DUT O ’ 
oy 
or 
cosw ot = Said 
a a -4 
On Sinw T F 4572) 


D 


555 


Again frequency warping 1s evident; but, as with the bilinear 
z-transform, prewarping may be used to produce useful results. 
>. A Summary of/S—Plane. to Z-Plane wanes tormactons 

Each of the transformations discussed in this section 
has advantages and disadvantages. In the design of digital 
filters the choice of a transformation must be dictated by 
these factors. 

First the standard z-transform, which is an expo- 
nential transform, requires the partial fraction expansion 
OE the transter functren H{s) and the attendant factouuZamuan 
of the denominator polynomial. The transform does preserve 
the shape of the impulse time response; but, because of 
frequency folding, useful results can be obtained only for 
bandlimited functions of the low-pass and bandpass varie 

The bilinear z-transform, being an algebraic trans- 
form, 1S easy to implement and requires no reduction of the 
transfer function ins. However, because of the frequency 
warping, this transformation should be used only for func- 
tions with piecewlise-constant magnitude characteristics. The 
bilinear form does preserve flat gain characteristics. 

The matched z-transform, like the standard z- 
transform, 1S exponential in form, but requires only the 
factored form of H(s), not the partial fraction expansion. 

It preserves the shape of the frequency response character- 
istics, and gives good results for high-pass and bandstop 


filters. 
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TLE. DIGWALA TE hRYDEOLGN TECHNIGUES 


It has been shown that digital filters may be described 


byacdigital transfer functions of—prhemtomn 





N -i 
5 a, Z 
= G2)... ee 
Hg(ez } is X (Zz) = Ph ie — 1) 
1 + ¥ b.z 
i1=1 


Equation (3-1) is perfectly general and applies to all 
Gigital-filter functions. The nature of the filter function 
is defined completely by the a. and ba, and the specifica- 
elLems of seees coefficients becomes almost the entire filter- 
design problem. | 

Digital filters are generally classified as being "recur- 
Sive" or “non-recursive." Referring to Eq. (3-1) if at least 
one of the b, are non-zero, the filter is said to be of the 
recursive type; 1.e., the computation of the present value of 
the output requires not only the present and past values of 
Gee input, but also the past values of the output (feedback) . 
If all the be are zero, then the filter is said to be of the 
non-recursive or transversal type; i.e., the output is simply 
a linear combination of weighted input values. This filter is 
also referred to as a "tapped delay line" filter. 

The design methods for the two classes of filters differ 
markedly and each class has its distinct properties. The 
non-recursive filter generally requires a large number of 
terms to obtain relatively sharp transitions in frequency 


response. However, the non-recursive filters have excellent 


S15 


phase characteristics, being quite linear 1n some cases. 
In contrast, the» recursiwe filter can provide sharp cutonrst 
characteristics with relatively few terms, but the phase 
characteristics are poor. 

In Section IV it will be shown that real-time implemen- 
tation of digital filters is presently much easier and more 
practical when the recursive filter structure 1s used. For 
this reason the discussion of non-recursive filters in this 
thesis is brief. All digital-filter design examples are 
recursive type filters to permit subsequent real-time 


implementation. 


AS ENON-RECURSIVE DIGITAL FILTERS 

The design methods for non-recursive filters fall into 
three categories. The first category 1s the classical 
numerical approach. In this method classical interpolation 
and differentiating formulas for equally spaced data can be 
digitalized directly by making correspondence between the 
unit advance operator "E" in numerical analysis and the 
z-transform operator. These filters are good only at low 
frequencies and they correspond to using the first few terms 
of a power-series expansion. This method is restricted pri- 
marily to filters whose function 1s integration, interpola- 
Chom, Gltrreremtlation, etc. 

The second design category involves the use of the Fourier 
series. If the magnitude of the desired frequency response 
is expanded in a Fourier series over Jw|<w/2 and the series 


i wWerelren in cermmaror cosines (or Sines) the result is 
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H(w) = 2 a, cos(nuT). (SicaD ) 
n=O 


But, since zt - a JUT, Eq. (3-2) can be written as 


a Ca eZ . (3-3) 


O 


lm 8 


uleze) = a + 5 
ig vl 


which is the required filter. The filter has an infinite 
number of terms and is realizable only if an approximation is 
made by truncating the series. The approximation is good 

for a rapidly converging series, but will be poor if there 
are discontinuities or steep transitions in the desired 
frequency response. 

The third category of design methods for non-recursive 
filters involve attempts to achieve a better convergent 
series, or to systematically weight error functions. Methods 
exist which utilize the Chebyshev series, least-squares 
approximation and "window" functions. Kaiser [1] treats this 
subject extensively and he presents several examples illus- 


trating the design methods. 


M. RECURSIVE DIGITAL FILTERS 
, 7 Nene een ee Pr ememmaneres tm 

Recursive digital-filter designs may be accomplished 
either in the time domain or in the frequency domain. His- 
torically designers of sampled-data filters were concerned 
with step-response, peak over-shoot, etc., and most transfer 
functions were analyzed and synthesized to satisfy time- 
domain criteria such as these. More recently, since the 


scope of digital signal processors has widened, more emphasis 


has been placed on steady-state response in the frequency 


5 / 








domain. Goldand—Rader-{4} have~formal:zed—various—proce - 
dures—for..digital-filter design in the frequency domain... 
Golden—([6]--hasalso..contributed to the.development of these 
metheds-;--furthermore-; he-has—written..a computer program 
whitch-essentratbiy~performs~-add..the..design computations to 
be-found™in=this~sectron. The-program—-was—not’ dvVarrabte—to 
the-author~of~thrs* thesrs+ 

There oe general approaches to the design of recur- 
sive filters. They are the direct deSign method and the | 
indirect design method. This latter method, which appears 
to be more popular, uses many of the standard continuous- 
time filters as intermediate designs. 

Pei ecr Design Method 

This method employs computational techniques to 

determine the coefficients, ay and bo, directly from the 
filter specifications. Some methods involve the use of 
polynomial-approximation theory. Another technique is to 
locate poles and zeros directly in the z-plane ee, and then 
possibly to perform an iteration or trial-and-error proce- 
dure. These design methods might be employed where the 
specifled response is not of a standard form (low-pass, 
bandpass, etc.). They—witt—net-be-dtrscussed further-tmn-thrs 
theese. 

2 Wileeated Fi eer Design Vraytialog-mMi leer" Design 

The indirect method of digital-filter design is 

accomplished in two steps. First an analog filter H(s) is 


designed to satisfy the filter specifications. Then a 


5.5 


transformation of the analog function is performed to déter- 
mine a digital transfer function H(z). In this manner use 

can be made of the many design and synthesis techniques for 
continuous filters which have been developed over the years. 


benwoi_trhe_extensive—-treatmert of thts 





subj ect_in_the..literature: 
a. Topics in Continuous-Filter Design 

Standard...forms.of-continuous—fitLters--have~been” 
developed~and~arewetldocumented’rnthe literature --s0-. 
The filters are referred to by the names Butterworth, 
eiebyshev, Elliptic, Bessel, Lerner and others. Each type 
of filter has certain characteristics which determine the 
most appropriate use for that filter. Of these many filter 
types the Butterworth and the Chebyshev designs are now 
briefly outlined to provide a basis for the digital-filter 
design examples which follow. 

(1) The Butterworth Filter. The Butterworth 


filter can be specified by the relationship 


jfoe (4 4) (9-4) 


1+ (w/w.)*™ 
where W, is the cutoff frequency. | H(jw) |? is the squared 
magnitude of the filter transfer function. It can be shown 
that the transfer function H(s) iSewa BAatLOnaIeiunctton of s 
with a constant numerator and a denominator polynomial of 
degree ns For example, a’ fourth-order Butterworth rilter 
with We normalized to 1 rad/sec has a transfer function 
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Figure 3-l. Response of a Butterworth Low-Pass Filter. 
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Figure 3-2. Response of a Chebyshev Low-Pass Filter. 
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LE | 
ey (3-3) 
s + 2.613s + 3.414s° + 2.613s.+ 1. 


The coefficients for other-order filters are readily found 
in tables GM#M™. Plots of the Butterworth frequency response 
for three values of n are shown in Fig. Arh. Generally the 
response is characterized by a flat passband and a stopband 
attenuation of 20 ndb/decade. 

(2) The Chebyshev Filter. The Chebyshev 


filter is specified by | 
Cl ~2) 
lH (jw) |7 = 5+; ’ (se) 
l+e V5 (0/4) 
where We is a Chebyshev polynomial of order n. The Chebyshev 
filter is characterized by an equal ripple in the passband 


and a monotonic decay in the stopband. The ripple amplitude, 


6, is given by , 
( 4-l 
6 = 1] - ————. 3-7) 


The transfer function of a Chebyshev 
falter is also a rational function of s. The coefficients of 
the transfer function are determined by the specified value 
of € and are calculated by a somewhat involved but well 
established procedure LO}. A plot of the response of a 
typical third-order Chebyshev filter is shown in Fig. 4-2. 
Generally the Chebyshev filter will have a steeper slope in 
the stopband than a Butterworth filter of the same order. 

(3) Frequency Scaling. Since most con- 


tinuous filters are designed using normalized frequencies 


Al 


mee. , WL = 1 rad/sec), frequency scaling iS required to 
shift the response of the filter up to the frequency band 


of interests» If we denote the desired cutoff frequency 


by w., then to shift the frequency response of the normalized 
transfer function to We make the substitution 
[5 
S n = om ° ( 3~6-) 


e 
(4). Frequency al sbanstormaeions. Higitpass, 
bandpass, and bandstop filters may all be obtained from a 
low-pass design by effecting the following frequency 


transformations: 


Low-pass-to-high-pass (1.60) 
k; | 
sot = (3=9e) 
Low-pass-to-bandpass 
Ww Ww a ) 
O S O | 


Low-pass-to-bandstop 


so —t_ (4-66) (3-9e) 


—- tye ES 
In the expressions BW = Oe Deal an i ¥ 2 ate Neel 
(the geometric mean). 

These transformations are all performed on 
the normalized low-pass H(s) since they contain their own 
scaling. Note that when using the bandpass and bandstop 
transformations Ww, 1s not necessarily in the center of the 
band. Consequently final designs may be slightly asymmetric 


relative to the intended center frequency. 
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b. The Design Procedure 

A desired frequency response 1s normally spec- 
ified by the nature of the characteristic in the passband 
and the amount of attenuation in the stopband. For example, 
it may be required to design a bandpass filter which is flat 
in the passband centered at 1000 Hz, down 3 db at 800 Hz 
and 1200 Hz (cutoff frequencies), and down at least 20 db at 
500 Hz and 1500 Hz. The flat characteristic would normally 
call for a Butterworth design, the 20-db attenuation spec- 
ification would determine the order of the filter and the 
3-db frequencies would specify the required frequency trans- 
formation. In the case where a final digital-filter design 
is required, the sampling frequency would be an additional 
and extremely important specification. 

The next decision to be made is the type of z- 
transformation to be used. Generally the standard z- 
transform would be used only for bandlimited functions to 
prevent frequency folding as discussed earlier. Bandlimited 
functions are typically functions of s whose denominator is 
of high degree and whose cutoff frequency is a small fraction 
Of the half-sampling frequency. Usually the success of a 
standard z-transform design cannot be assured without 
actually testing the final filter. 

The bilinear z-transform works well when the 
frequency characteristic is relatively piecewise-constant. 
Hence, most Butterworth designs and small-ripple Chebyshev 


designs could be transformed Satisfactorily. The algebraic 
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form of the bilinear z-transform makes it easy to implement; 
but prewarping of the critical frequencies in the s-domain 
is required to achieve the desired response in the z-domain. 
Finally the matched z-transform may be used to 
achieve a frequency characteristic that is high-pass, band- 
stop, or other than piecewise-constant. 
The digital-filter design procedures discussed 


above will be illustrated by examples in the next section. 


Pe EXAMPLES OF RECURSIVE DIGITAL FILTERS 

The following examples were chosen to illustrate the 
digital-filter design methods which have been discussed. 
The order of the filters was kept low to facilitate the 
hanging Of the numbers. All calculations were persorned 
on an electronic calculator, and the results were rounded 
to the sixth decimal place. Fex-more-cemplex—filters the 
use _of a—ecomputes—_program-such—as—the. one-wertten—by=Golden 
[6)}>would--almost—certainly be reguixed. 
EXAMPLE 1. Design a low-pass digital filter to operate ata 
Sampling frequency of 2000 Hz. The cut off (3-db) frequency 
is to be 200 Hz and the response is to be essentially flat in 
the passband. 
1. A second-order Butterworth filter is chosen with a sub- 
sequent transformation using the bilinear z-transform. 


Prewarping of the analog cutoff frequency 1s required. From 
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2D AO eae Cae yer 3914159) = 0.324920. 


—— ~ 


Note that the analog cutoff frequency Wn has been prewarped 


to (0.324920) (2/T) = 207 Hz. 


2. A normalized second-order Butterworth filter is MA) 
oe | a aaah 
Siwtteed M92 de3seck “L ' 
Substituting’ s = s/0.324920 igng@yPtSrrr gives 
_™ 
H\(s) = : 0 2105 Si 3 ; 
Ss“ + 0.459506s + 0.105573 LEOY 


1 
3. The bilinear z-transform is taken by substituting 


s = (z-1)/(z2+l) in Eq. <@™. The result is 


W. 
H(z) q 0.105573 (2° +22+1) , 


2° -22+140.459506(z=1) (+1) + 0.105573 “e-Hoeea) 
which reduces to 

oemewass Hom 34e10 2 2 + oMe7455 2 - 
H(z) = se teee 7 ees eGe.-3,) 


1 = 1.142980 z ~~ + 0.412801 zg 
The frequency response of this digital filter was evaluated 
using Computer Program 1, and the results are plotted im abibe & 
Poyseos>. The 3-db frequency is exactly 200 Hz. The poles 
o Fire progrom 
# The factor 2/T is dropped here and in the subsequent 
Substitution for s to avoid handling large numbers. Should 


the factor be retained the result can always be reduced to 
the same final answer. 
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of the transfer function, @@~- (GMB) are: 
Zz = 0.571M90f 702293598 
These poles are well inside the unit circle indicating a 


stable filter. UIT aoaceeeanietainnareseenettlinasinane chai ee 
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EXAMPLE 2. Design a bandpass digital filter to Operate™at™a 
sampling frequency of 2000 Hz. The cutoff frequencies, Wry 


and w are 200 Hz and 500 Hz respectively, and the response 


Dy 
is to be flat in the passband. 

1. A second-order Butterworth iS again chosen to be used 
with the bandpass bilinear z-transform given in Eq. (2-41). 
Mais transformation will yield a fourth-order functeon in 

z for a second-order function in s. 


2. Equation (2-42) is written again for convenience. 


cosw T - cosw_T 
oO D 


me ¥ SinwpT (3-14) 


It is desirable to have the analog frequencies associated 
with each of the critical digital frequencies be negatives 
of each other. Then the digital transfer function which 


results by transforming H(s/|w,|) will satisfy the response 


specifications at both critical frequencies, On and Woo ° 
To accomplish this solve for cosw T in 
cosw it = cosuy,t » = cosut ta cosuy ae, 
- Sinw,,1 S1nw,5T 


After considerable manipulation of Eq. (3-15) using trigo- 


memetric identities, thaeresult iis 
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lee | 
COS. = Uae Fite Ie 
cosa. T.= 4 we Be ” (3-16) 
a a (u=ar — O-nT” 
ee p2-’ 
Equation (3-16) need be derived only once. It then becomes 


a formula for use in the design procedure for the bandpass 
bilinear z-transform. 


3. Returning to the example 
cos & (700 - Wn - 7} wa 
cosw oT = 2 EEE 0.509522 a Cei7 ) 


cos = (-300 - 2n « T) 


4. Then, using Eq. (3-14), 


Woz T = (200) (2m) (0.5.x 07") = qe cage, 
and 
_ 05509503, =acosstmerss2 : 
4 = =  -£rar O1Gpeeo ie) OSS 2 6 (3 18) 


An identical result is obtained for “..T = 1.57080 verifying 


D2 
tite valadity of Eq. (3-16). 


5. Substituting s = s/0.509523 into the normalized second- 


order Butterworth transfer function gives 


= 


H{s) = 


0.259614 (3-19) 


B° F Ob 20e746 + 0.259614” - 


a 


6. The digital transfer function is now obtained by substi- 


Z 


ether aicmem——ize — La019046z = eee =e) j)cilteomnae (3-19). 


The result after considerable algebra is 
(3-20a) 
0.131106-0.8622122 -+0.1311062 — 


ar a ee 
1-1.400064z ~+1.272216z “-0.658419z ~+0.2722172 


oT 
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9 SS wn (3-20b) 


(0 .362085-0.7241702 +0 .3620852z ~)° 
2 


(1-1.222023z "+0 .603825z “) (1-0.17804127++0.45082127 ) 


\ er 
Nes § x, \ &. ~ = 4 - 


Hi{z)= 


The frequency response of the filter is evaluated using Eq. 
(3-20b) and the computer program. The results satisfy the 
specifications exactly and are plotted in Fig. 3-4. The 
poles of this bandpass digital filter are located at 


0.611011 + 40.480094 
0.089020 + 40.665505 


Z 


Z 
Again the filter is stable. 
EXAMPLE 3. Design a bandpass digital filter with a frequency 
characteristic within the passband equal to the reciprocal of 
a gaussian or normal curve. (An application of this filter 
is given in Section VI, C.) Refer to Fig. 3=5.°' The center 
frequency is 450 Hz and the sampling rate is 2000 Hz. 
1. The shape of Fig. 3-5b may be approximated by a second- 
order Chebyshev filter. (An nth-order Chebyshev filter has 
n peaks and n-l troughs. Hence, a second-order filter with a 
large € will have a dip in the center of the passband.) The 
design of this particular Chebyshev filter will only be 
summarized here. 

A value of € was calculated to give the dip a shape which 
approximates the reciprocal of a normal curve. To do this a 
relation between the statistics of the curve and frequency 


had to be aaa. ened.” The normalized filter was then 


: In Section VI, C, the origin of the gausSian curve is 
seen to be the spectrum of a CW carrier which is frequency 
modulated by noise. Hence, the relation between the spread 
of the curve and frequency is determined by the actual 
hardware being used. 
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a. Gaussian Distribution in Frequency Domain 
yw = mean = center frequency 
Oo = standard deviation 


-20 -O u O 20 


b. Reciprocal of Gaussian Distribution 


fs 
O 
| ritter Bandpass 


fol fFo2 


Figure 3-5. A Noise Filter Specification. 


Bye 


designed using standard methods [10]. Final'y a frequency 
transformation was performed to give the required bandpass 


filter, =the result ss 


H(s) = 0.515422s7/(0.000000052s* + 0.000039934s° 

5 (3-21) 

+ 1. 2501L15878s" +°2035 7773 Mibs + 26704094 G54).. 

The frequency response of this filter,°|H(jw)|,is shown in 
Fig. 3-6. The results are quite satisfactory. The problem 
at hand is then to obtain a digital transfer function which 
will yield the same frequency characteristic. 
2. The matched z-transform is chosen to obtain the required 


digital transfer function. The poles of the transfer func- 


Ebnem H (sis ane 


-288.74512 + 34639.3984 
- 95.23544 + 31521.2546. 


S 


S 
Digital transfer function coefficients can then be evaluated 
using the matched z-transform relationships given in Eq. 
(2-40). The results for the two second order transfer 


functions are 


ay = 1S 8. 

a, = -1.0 

by =@ 1.178165 
b. = 0.749203 

and 

ag = iP. 0 

a, = -1.0 

by = -1381435 
b. = OOO . 


6 Figure 3-6 is the actual output from the Calahan net- 
work analysis computer program [ll]. This program is an 
extremely useful aid in the present context of analog 
transfer function synthesis. 
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Then the complete digital transfer function is 
= _ (3=22) 
K(l=2 Hye z) 


fez) = al . 2 -—] -—?) ! 
(1+) 2173 065 Zee On 4920 37) 112.3 Bas ees 091.5977 ) 


where K is a gain factor. 

3. An examination of the frequency response of the transfer 
function given by Eq. (3-22) reveals a difference of 5.5 db 
in magnitude at the two peak frequencies. If a first-order 
zero at the half-sampling frequency 1s inserted in the 
transfer function, the desired symmetrical response can be 
better approximated. The response of this modified transfer 
function is plotted in Fig. 3-7 where the gain factor K was 


| adjusted for 6-db attenuation at center frequency. 


iS ya) 


IV. IMPLEMENTATION OF DIGITAL FILTERS 


Let H(z) define an arbitrary digital transfer function 
derived by one of the methods discussed in Section III. The 
solution forms thesoutpume tlansié Orange Y ( Zz) gee Bibz ea) Ganmlse 


Wigeeiiee Wl ees 


b. z > Y¥(z) ; (4-5 


or, in the time domain interpreting i as the unit delay 


operator, 
b.  -f (4-2) 


This linear difference equation can be realized directly 
using digital components. Refer to Fig. 4-1. Each rectangle 
inscribed with = represents a one-sample delay. In real- 
time hardware this means a one-word memory whose capacity 

in bits is determined by the number of bits of accuracy used 
to represent x, and y,. Also a multiplier is required for 
each branch inthe filter. In general these are variable- 
coefficient multipliers to permit changing of the transfer 
function. Notice that these multiplying coefficients are 
the same coefficients which appear in the transfer function 
H(z). This is the basis for the earlier assertion that the 
design of digital filters is essentially complete once the 
coefficients have been specified. This is in sharp contrast 


to analog-filter desiaqn where the determination of the 
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Figure 4-1. A Recursive Digital Filter. 


De), 


transfer function H(s) is often only a small part of the 
design problem. 

The adders and multipliers together comprise the arith- 
metic unit of a special-purpose computer or digital filter. 
The implementation of Eq. (4-2) would involve clocking this 
computer at a rate f=1/T, and performing the indicated 
Calculations during each sample interval. After 
the calculations are complete, the contents of the storage 
registers, or delay elements, must be shifted one location 
to be in the proper position for the next sampling period. 
For example, if the input value at time t = nT is stored 
as XL, then this value must be stored in the location for 
X 1 One period later at t = “oe ae 

In general, the expression for H(z) = Y(z)/X(z) is a 
ratio of polynomials in 2 where the numerator is of degree 
N and the denominator is of degree M. If N > M then the 
computation of the current value of the output x, in Eq. 


(4-2) would require knowledge of the input value x For 


ii 
example, if N=5 and M=4, then the computation of ya would re- 


quire the value of x, among others. Since this is physically 


5 
impossible, it can be stated generally that for a digital 
transfer function to be physically realizable the degree of 
the numerator must be less than or equal to the degree of 
the denominator. 

It is seen that the realization of digital filters in- 


volves basically a computational device to solve the linear 


difference equation which is defined by the digital transfer 
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function of the filter. However, the sequence of performing 
the arithmetic operations and the method of storing inter- 
mediate answers do have an effect on the efficiency of the 
filter and the accuracy of the results obtained.’ A digital 
meecer 1s said to be in canonical form if, for a given order 
n, a minimum number of adders, multipliers and delays are 
utilized. Three canonical filter structures are now 


discussed. 


A. THE DLRECT Bon” 
Consider the general second-order digital transfer 


fnnctLon 


H(z) = —————_—__———_—_——_ .. (4-3) 
Ll + by 


H(z) can be written as the product of two transfer functions, 


H, (2) and H,(Z) . Then Eq. (4-3) becomes 
aia. a ee 
cz) 1 eo: a ee Se 
X(z) L+b, 2 b+ pb, 2°? 
2 
and 
wz) = X(2)H, (Zz) H(z) : a) 

Define the intermediate variable W(z) = X(z)H, (Zz). Then, 


ar the usual notation 


7 Efficiency refers to the speed of computation and the 
effective utilization of hardware. Both these items have 
become important "yardsticks" for measuring the usefulness 
@® a digital fxrlter- 
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ig n 1 n—-i. Z 
(a6) 


= a WwW + a.w + a 
Yn 0) Sal di 


The general result is that the single difference equation of 
the form of Eq. (4-2) is converted into two difference 


equations of the form 


M 
Wy Se EP ies 
T=. 
Ca.) 
N 
ve = boa ee ae 
iL=6, 
The schematic representation of Eq. (4-7) is shown in Fig. 


4-2. Notice that Fig. 4-1 and Fig. 4-2 are similar in form 
but Fig. 4-1 contains redundant delay elements. Hence, Fig. 
4-2 1S a canonical filter structure. This means a more 
efficient use of hardware. 

Any nth-order digital filter may be implemented in the 
forms of Figs. 4-1 and 4-2. These are both called direct 
forms and are the simplest of the recursive filter structures. 
However, because of numerical considerations, the direct forms 
are normally satisfactory only for low-order filters. These 
numerical considerations will be discussed in more detail 


mamoection IV, Ds 


B. THE CASCADE FORM 
The general nth-order digital transfer function may be 
factored into a product of second-order transfer functions. 


Mie result is 
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=i -2 
(nn os 2 O52 + Oo. Zi 
H(z) =a, 1 ——-—.———___, , (4-8) 


O- | = 1h 
i=a 1 ft Bay Z. + Bos Z 


where™m is them@integer part of @(n+P)/72. Netemch aceror’n (27 
with poles and zeros in complex-conjugate pairs the coeffi- 
cients a and 8, will always be real. The resulting 


difference equations are of the form 


Y, (2) = H, (2) X(2) 
Yo (2) i H, (2) ¥, (2) 4 e 
Y (zy - Het 2) Y 2 a 


by the same arguments which were made in the preceding section. 

The schematic representation of Eq. (4-9) 1s shown in Fig. 
4-4 It 1s seen to be a cascade of second-order subfilters 
and this structure has several important features. First 
the structure can easily be multiplexed. By uSing a basic 
second-order filter structure and appropriate timing and 
control circuitry, a high-order filter can be implemented by 
time-sharing the second-order structure. The appropriate 
coefficients, Os and Ba, are supplied for the m subfilters 
in time sequence. 

Also notice that if H(z) contains Zeros on the unit 


Circle, sthescockLicrent > in the subm@ilter which contains 


PA 
a second-order zero will be unity. This means that one less 
MULE edeieC Aili Dele WeCc due. pxGliGbed eodiwatLon OGGgues Lige— 


quently. when the bilinear z-transform is used to design band- 


pass and bandstop filters. 
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H, (z) H(z) ae H (2) Yn 


Figure 4-3. The Cascade Form 








Yo 
co 
x —- Yn 


age) 


Figure 4-4. The Parallel Form 


Gye LHe, PARALLEL FORM 
The third canonical form*®et a, ,digwea falter is obtarmea 
by expanding H(z) in partial fractions of second-order terms. 


The result.1Ls 


-i 
m \ius cee Ye Z 
| > Cea (4-10) 
i=l l+ oes -y Pontes 


where ‘ee 2 yh am Equation (4-10) yields difference equa- 


tions of the form 


Y, (2) ~ ¥o% (2) 
Y, (2) = Hy, (2)X(z) 
: (4-11) 
ea) = HH 62) X(Z) 
m 
Y(@ - ee LZ) ; 
i=0 


and the filter is sgructuré@ as snown in Fig. 4-4. The 
parallel form is comparable to the cascade form in perform- 
ance. The choice between the two might be made on the basis 
of the amount of algebra involved. Depending on the z- 
transformation employed the amount of work in arriving ata 
final filter structure can be considerable, especially for 


high-order filters. 


8aiternatively the expansion_may be made of H(z) written 
as a Eunction- of Z ratmer than z 
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D. NUMERICAL CONSIDERATIONS 

Until this point it has been assumed that both the cal- 
culation of digital filter coefficients and»the manipulation 
of these coefficients in the filters themselves were done 
with infinite accuracy. Since this cannot be the case, and 
Since even the input samples to the digital filter are only 
an approximation to the instantaneous values of the analog 
Signal, the actual results may differ markedly from the 
theoretical results. 

As with any numerical process the accuracy of the cal- 
culated "answers" using digital filters 1S a very critical 
function of the hardware used and the computational technique 
employed. The finite computer word length is the basic 
limitation on accuracy. Errors introduced because of this 
limitation are: (1) finite precision in®filte® Goefivorents; 
(2) quantization of input data to a specific number of bits; 
aaa (3) xound-off error in multiplications and addations® 
References 1 and 2 contain rigorous treatments of quantiza- 
tion effects in digital filters. The results are stent 
below. 

Numerical problems can begin even with the digital filter 
design. The design methods of Section III involve factoriza- 
tion and partial fraction expansions. Generally it is 
recommended to perform these operations in the s domain and 
to digitalize the resulting subfilters. Although factoring 
of the polynemials ins Ws not required fer the biaiinear z-— 


taansform, better mesults arewobta@ned if thas is done. * ie 
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generation and manipulation of high-order polynomials in z 
sheuld be avoided. Inepartmemlar, finding roots of poly- 
nomials im z can be quite difficult because of the tight 
clustering of poles and zeros in the z-plane. 

The popular notion that a higher sampling frequency 
always produces better results is not true. Actually a higher 
sampling rate requires greater accuracy of computation. Per- 
formance may deteriorate if the sampling rate is increased 
beyond certain limits. 

Often the most critical numerical problem in digital 
filters is the quantization of the filter coefficients. 

These coefficients determine the locations of the poles and 
zeros for the filter and any perturbation in the coefficients 
causes a shift in the pole and zero locations. for high- 
order filters the poles are usually clustered ina small area 
near the unit circle in the z-plane. Perturbations in these 
coefficients may actually cause the filter to become unstable. 
Furthermore, the pole locations are much more sensitive to 
coefficient perturbations when using the direct filter 
structure. Consequently high-order filters are almost always 
realized in the cascade or parallel form. Weaver et al [7] 
present a procedure for determining the bounds on the sensi- 
tivity of pole positions to coefficient perturbation. Given 
these bounds the required coefficient accuracy can be cal- 
culated. 

The quantization of the input signal in the analog-to- 


digital converter is treated as a noise input to the filter. 
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The round-off error in the arithmetic unit is also treated 

as noise but its ehieryemtoethe £1 ltenemayeoecumminewanweus 
places depending on the realization structure. In the direct 
and parallel foriis’the"quantization noise due, to roundsoff 
error is simply additive. However in the cascade form round- 
off error may cause unsatisfactory operation because of the 
multiplicative effect in a serial system. A complete anal- 
ysis of various implementation schemes can be accomplished 


using noise models such as those developed by Gold and Pader 


[4]. 


E. HARDWARE CONSIDERATIONS 

Any device which operates on signals in real-time accord- 
ing to a digital transfer function is, in essence, a special- 
purpose aigital computer. The computer program, or 
equivalently the transfer function and implementation scheme, 
can be hardwired directly into the computer. The input data 
would then be the signal samples and the filter coefficients. 
It is anticipated that LSI will provide the mechanism for 
incorporating such a computer into a package that is small 
in size, economical in cost, and reliable in performance. 

The design of such a package is a formidable task, and several 
of the more important considerations are now discussed. 

The first consideration is the choice of the sampling or 
clock frequency. Unfortunately the requirements of a digital 
filter in this respect are self-opposing. To obtain sharp 
cutoffs, more narrow bandwidths, deeper notches, etc., 


higher-order filter functions are required with 
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correspondingly large computation times. ‘Since the computa- 
tion must be accomplished in a single sampling interval, the 
sampling frequency must go down to provide a longer time T. 
However, as the sampling frequency is lowered, the usable 
frequency spectrum is reduced because the relation W <W/2 
must be maintained. 

The second major conSideration 1S computer word length. 
The cost and complexity of a real-time digital filter varies 
in direct proportion to the number of bits which must be 
stored and manipulated. It appears that the key to success 
in reducing computer word length lies in developing new 
implementation schemes which are less sensitive to numerical 
inaccuracies. 

In general, therefore, 1t 1s of the utmost importance to 
design and organize special-purpose hardware which can per- 
form digital-filter type computations efficiently, accurate-- 
ly, and as fast as possible. One of the more demanding 
computations in a digital filter is multiplication. In 
order to illustrate some of the problems in logical design 
peculiar to digital filters a binary Serial multiplier Ws 
outlined in detail in Appendix A. The multiplier fdllows the 


general design proposed by Jackson et al in Ref. 8. 


F. THE FAST FOURLER TRANSFORM 

The methods discussed thus far in implementing digital 
filters all involve the technique of interpreting zt as the 
unit delay operator and obtaining a solution of linear 


difference equations directly in the time domain. This 
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technique ws) pertecrly wale fer vanally zing thersite aay—state 
response of anette transfer functions to sinusoidal inputs, 
and has been shown to be very practical in real-time appli- 
canienee However, for complex filters, and in particular 
for non-recursive filters with a large number of terms, the 
number of multiplications and additions tiay’ be préhibitive. 
Mueechese cases “filter implementationm™py higuespeee Conve 
lution using the fast Fourier transform (FFT) becomes 
desirable. 

Recall that the discrete convolution of an input sequence 


with a filter impulse response is given by 


~~ Swen x, (4212) 


where M is the number of samples being convolved. When M 
is large evaluation of Eq. (4-12) by z-transform methods 
becomes difficult and time consuming. 

Alternatively the sequence y, may be evaluated by (1) 
taking the Discrete Fourier Transform (DFT) of the input 
sequence and the impulse response of the filter, (2) multi- 
plying these two transforms together, and (3) taking the 
inverse DFT of the result. This procedure is very similar 
to uSing the conventional Fourier transform to analyze 
continuous-time signals in the frequency domain. However 


the theory of the DFT is somewhat more involved. 





| 9 This is due to the fact that the product of two DETR s 
in the frequency domain is equivalent to circular convolu- 
tion in time, not linear convolution as intended in Eq. 
(4-12). Nevertheless, the general results as stated above 
ave COrrece. 
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Thesepractircalptys.ofethasemethod laiessinathe abi lityeto” * 
use the FFT techniques to evaluate the DFT's. The fast 
Fourier transform 1s described in detail in Refs. 4 and 
12. Briefly it makes use of the fact that under certain 
conditions the calculation of the coefficients can be done 
iteratively with a considerable saving in time. Generally 
the value M (the number of sample SASS must be a highly 
composite number (contains many factors). When M is a power 
of 2 the FFT requires a number of computations proportional 
to Mlog.M vice mM? for the conventional convolution. 

Presently a large general-purpose computer 1s required 
to perform high-speed convolution via the FFT. It is antic- 
ipated that this algorithm will eventually be programmed into 


Special-purpose hardware for real-time applications. 
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V. REAL-TIME SIMULATION OF DIGITAL FILTERS 


In recent years large general-purpose computers have been 
used extensively to perform the computations required in 
digital eau processors. Computer Program 1 is an example 
of the frequency-domain analysis of a digital transfer 
function. In the time domain the calculation of an output 
sequence of a digital filter can also be programmed for a 
general-purpose computer. In this case the actual difference 
equations are solved and the program becomes a simulation of 
the digital filter. One limitation in such a scheme is 
immediately evident; a large computer memory is required to 
store the input and output number sequences. Nevertheless, 
such computer simulations of digital filters do find wide 
application. Indeed, most algorithms involving the use of 
the fast Fourier transform are performed on a general-purpose 
computer. 

With the increased emphasis on the real-time implementa- 
tion of digital filters it becomes desirable to be able to 
Simulate digital filters in real time. Such a simulation 
would provide an experimental basis for establishing specifi- 
cations for the realization of digital filters in special- 
purpose hardware. It would provide a means for investigating 
and verifying new design and implementation techniques. And 
finally it would provide a real signal to observe, measure, 


and analyze both quantitatively and qualitatively. 
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A real-time® simulationw~of “dtogital *f11lters must contain 
certain features 1f it is to be used as an experimental anal- 
ysis and design tool. First it must be possible to control 
the two basic design parameters--clock frequency and computer 
word length. Precise control of the sequence of arithmetic 
operations is required to simulate different filter struc- 
tures. And finally the simulation scheme must contain the 
necessary input-output interfaces to permit the complete 
processing of a signal waveform in real time. In short, it 
is required that a simulation scheme be developed which can 
duplicate as closely as possible an arbitrarily designed 
special-purpose computer or digital filter. 

A real-time digital-filter simulation in the above context 
cannot be accomplished using only a large general-purpose 
computer. Aside from the obvious lack of any analog signal 
interface, the control of computer word length and precise 
arithmetic sequence is difficult, if not impossible, to 
achieve in most large-computer operational environments. 
Furthermore, the use of a large computer for relatively simple 
computational algorithms for an extended period of time is 
extremely inefficient and, hence, undesirable. | 

On the other hand, the real-time simulation problem is well 
Suited to a scientifically oriented hybrid-computer system. 
Some of the characteristic features of such a system are: 

1. Built-in interface between analog and digital computers 


including analog-to-digital converters and digital- oe 
analog converters. 


2. Readily accessible analog functions such as amplifica- 
Eionpescatame, Clipolng, mmteqratren, ete. 
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3. Readily accessible digital functions such as clocking, 
COUNTING Scoring, Ca lel lating , ois@c calm ine), selec . 


4. Accessibility of such peripheral equipment as oscil- 
loscopes, vacuum-tube voltmeters, signal and function 
generators, line printers, card readers, tape recorders, 
and teletypewriters. 


5. Typical "hands on" operation of the computer system 
vice the "batch processing" type operation of a computer 
facility. 


In the following sections a real-time simulation for dig- 
ital filters is presented. The specific simulation scheme 
was designed for use on the hybrid-computer system located in 
the Department of Electrical Engineering at the Naval Post- 
graduate School. While the ensuing discussion must be 
specific in certain areas, the solution to the overall prob- 


lem 1S quite general. 


A. THE HYBRID-COMPUTER CONFIGURATION 


1. The Physical System 
The hybrid-computer system at the Naval Postgraduate 
School consists of the following major subsystems: 


1. Xerox Data Systems XDS 9300 Computer 19 
2. COMCOR Ci-5000 Analog Computer 


3. Hybrid interface subsystem consisting of: 
a. Adage VMX multiplexer 
b. Adage VT-13 14-bit A/D converter 


some of the characteristics of the XDS 9300 digital 


computer which are of specific interest are: 


ley 24-bit wordwincitidimg sitgm birt 
2. META-SYMBOL Symbolic Assembler 


10 This nomenclature was formerly "Scientific Data Sys- 
tems SDS 9300 Computer." 
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1.75-microsecond cycle time 

1.75-microsecond fixed-point add time 
7.0-microsecond fixed-point multiply time 
direct memory access for A/D converter output 
integral 15-bit D/A converter 


3 index registers available to the programmer 


o OI OD MO F&F W 


priority interrupt system 

The Ci-5000 analog computer has a full compliment of 
operational amplifiers, potentiometers, passive analog com- 
ponents, logic components, clock sources, and computer 
control devices. These items will be explained as necessary 
in the discussion. 

The Adage A/D converter produces a 14-bit binary 
fraction including sign bit from an analog input signal with- 
in the range + 100 volts. This 14-bit fraction is stored Tm 
the most significant (leftmost) portion of an XDS 9300 24-bit 


memory location. That is, 


sO ON 3777/6000 


8 
Om. = 00000000, 
SOOTY, = 40000000, 


where the decimal point is assumed after the sign bit, and 
negative values are stored in two's complement form. 
2. The General Simulation Scheme 
The digital-filter simulation which follows is built 
around a digital-computer program written in a symbolic 
programming language. This symbolic language, called META- 
SYMBOL, is the means whereby a machine-language program is 


generated in the XDS 9300 computer, 
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A machine-language program is desirable for several 
reasons. First it provides the rigid computer control and 
aecessibility! that is required for real-time filter simu- 
lation. For example, arithmetic operations such as addition 
and multiplication are specified precisely in the desired 
sequence. The location in memory of each variable, machine 
instruction and data word is known, and the flow of digital 
information in the computer is completely specified by the 
program. Furthermore, a machine-language program permits 
the operator to monitor the execution, enter the program to 
insert modifications, or to stop the program to examine 
intermediate results. In general, the programmer not only 
tells the computer what to do, but how, where, and when to do 
it. This complete control of the situation is required for 
the accurate simulation of a special-purpose computer. 

The analog signal which is to be filtered is fed to 
the patchboard of the analog computer. It is amplified, if 
necessary, and then fed directly to a trunk line that is 
connected to one of the input channels of the A/D converter. 
Operation of the A/D converter 1s under control of the dig- 
ital-computer program. At a certain point in the program the 
analog signal is sampled and the quantized signal magnitude 
is stored directly in a previously specified location in 
memory. The digital computer then begins to perform the 
arithmetic operations of the digital filter. The filter. 
coefficients had previously been stored in memory and the 


sequency of operations is programmed to follow the filter 


is 
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implementation scheme exactly. Upon completion of the cal- 
culations the output value for that sample period is stored 
and also made available to the D/A converter. When the D/A 
conversion is complete the filtered analog signal is fed back 
to the analog patchboard where it is rescaled and applied to 
an oscilloscope or spectrum analyzer. Back in the digital 
computer the various input and output sample values which 
are used for computation are shifted in memory to be in the 
proper location for the next sample period. The computer 
then waits for a clock pulse to begin the next iteration of 
the entire sequence. 

The master clock is located in the analog computer. Each 
clock pulse causes an interrupt in the digital computer 
Operation and immediately causes it to go to the starting 
point in the iterative program. Obviously the digital com- 
puter must have completed the calculations for a given sam- 
pling period before the next clock impulse occurs. Thus the 
computation time of the digital filter sets an upper limit 
on clock frequency. 

3. The Hybrid-Computer Diagram 

A block diagram of the simulation scheme is given in 
Fig. 5-1. An audio signal generator is used as a signal 
source. Because of the +l00V input range in the A/D con- 
verter the original signal is scaled by a factor of 10 in 
amplifier A034. This means the input signal must now be 
mem be0. tous 10.0v. The unscaled signal, Cans ise also” fed EO 


the upper trace on the oscilloscope. 
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The output of the D/A converter, which is the filtered 
signal, is fed to a unity-gain amplifier A036 through poten~ 
tiometer PO36. Assuming no additional scaling of e., was 
required in the digital computer a setting of 0.1 on P036 
should cancel the gain factor of 10 on the input side, and 


the magnitude of e should be correct according to the 


out 
transite rwiuietion,of thesfilter. If a given vakwe ef ern 
produces overflow in the computer's arithmetic unit, the 
ee im@icatom-onethesconsOleéeiselit.@ Inethas casesmtit= 
signal input level can simply be reduced, or additional 


z OG Ne can easily be accomplished in 


scaling (usually 2. 
the computer program. Potentiometer P036 is then adjusted 


to equate the net scale factor of the entire loop to unity. 


B. THE DIGITAL-FILTER SIMULATION PROGRAM 

Computer Program 2 is a general program which can be used 
for the real-time simulation of digital filters on the hy-~- 
brid-computer system discussed in the preceding section. 
The source program 1S written in the META-SYMBOL language 
and appears on the right side of the page. The block Cu 
numbers on the left side of the page is the corresponding 
machine-language program generated by the META-SYMBOL assem- 
bler in the XDS 9300 computer. A complete description of 
the operation of this computer and the mechanics of the META- 
SYMBOL language will not be given here. Detailed information 
on both subjects is contained in Refs. 13 and 14 respectively. 
However, the following limited information is provided in 


order to explain the simulation program. 
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1. Baé@kground™igtciiatiromomsthe {DS .930 0eC@omputer 
The META-SYMBOL source program is originally punched 


on cards. Refer to Computer Program 2. Card column 1 is 
printed in the Location of Sthe@asterisk¥im Dime numberela 
The format for the symbolic program is as follows: 


1. Card columns 1-7 contain the label field. This field 
is used for address identification or labeling. 


2. Card column 8-15 contain the operation field. Either 
mnemonic machine instructions or assembly directions 
are entered in this field. 


3. Card columns 16-30 contain the operand field and may 
be used for various purposes. Put very simply, the 
operand answers the questions "what?" or "where?" 
regarding the operation field. 


4. Card columns 31-80 contain the comments field and are 
ignored by the computer. 


The computer word length in the XDS 9300 is 24 bits 
Gr eight @etads daigats.. The smaghtmost.e1ghth digits in the 
block of numbers in the program printout represent the ma- 
chine instruction corresponding to the META-SYMBOL instruc- 
tion on the same line. Each machine instruction and data 
word in the program is stored in the memory location given 
by the first five digits in that line. For example, refer to 
line 34 in the program. The machine instruction 01600132 is 
stored in memory location 00035. Furthermore, the eight 
octal digits of the instruction word are formatted as follows: 


Digit 1 contains address and index register information. 
Digits 2 and 3 contain the instruction code. 


Digits 4 through 8 contain the address of the operand. 
An example will help clarify the discussion. Again referring 
to line number 34 in the program, the META-SYMBOL statement 


reads, "LDA «+** SUM+1l." This means simply, "load the A 
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register with the contents of memory locatiqn SUM+1." The 
operation code for LDA is 16. Therefore, the machine in- 
struction states, "load the A register with the contents of 


memory location 00132." From line number 100 it is seen 
that the location in memory which is one greater than the 
leocatweon | resexmved sfonseSUM=ers4nefactshocahmons number 00032. 
Address modification using indexing and indirect 
addressing is somewhat more involved and will not be ex- 
plained here. Moreover, many of the instructions in the 
program are required simply to make the program work and do 
not warrant explanation beyond the comments in the printout 
itself. Most of the "subroutines" (ADOPS, etc.) are re- 
quired by the hybrid interface software. 
2. Background Information on the META-SYMBOL Language 
The most important part of the simulation program 
involves the use of the computer's arithmetic unit to 
perform the prescribed calculations on certain sample values. 
The following list contains selected machine instructions 
from the META-SYMBOL program language which are used to 
Simulate digital filters. 
Pmseruction Translation 


LDA Load the A register with the contents 
of the effective memory location. 


STA Store the contents of the A register in 
the effective memory location. 


ADM Add the contents of the A register to 
the effective memory location. 
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MUL Multiply the contents of the A register 
by the contents of the effective memory 
location. The full product is contained 
in the A and B registers. 


ALSD The contents of the A and B register, 
treated as a double-length register, are 
shifted left a specified number of bits. 


ARSA A right shift similar to that above 
invelving the A register only. 

BRM Mark Place and branch. 

BRU Branch unconditionally. 

LDX Load the specified index register with 
a Specified number. 

BRX Branch only if the associated index 
register has not reached its terminal 
state. 


The general use of the above instructions in the digital 


Simulation are: 


i, 


Zs 


3. 


LDA, STA, MUL, AND ADM are used to add or multiply and 
to store results in a desired location. 


ALSD is equivalent to multiplying by 2” where n is 
the number of bits shifted. 


ARSA is equivalent to dividing by ee 


BRM is used to shift computer control to a desired 
"Subroutine." The return to the program is automatic 
when execution of the subroutine is complete. 


BRU is equivalent to a FORTRAN GO TO statement and 
is used in this program to iterate a sequence of 
instructions indefinitely. 


LDX and BRX provide a "DO LOOP" mechanism whereby a 
certain group of instructions is performed a pre- 
scribed number of times. Indirect addressing is 
often used to shift an effective address for each 
new iteration. 


Fixed-Point Arithmetic 


The XDS 9300 computer performs either floating-point 


arithmetic or integer arithmetic. Since floating-point 


arithmetic is far more complicated, and since it would almost 


never be used in special-purpose hardware, integer arithmetic 


is the most practical alternative. However, in the digital 
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filter itself, fixed-point arithmetic is desired whereby 

a mentally placed decimal point in a computer word remains 
fixed throughout the calculations. Hence, a conversion 
scheme 1S required. 

-Recall that the basic word length in the XDS 9300 
computer is 24 bits including sign. The output of the A/D 
converter is a 14-bit word corresponding to a binary frac- 
tion of magnitude less than one. That is, the first bit is 
a Sign bit and bits 2-14 are magnitude bits assuming a 
decimal point after the sign bit. The least significant 10 
whits Of a sample valuc in memory, therefore , amemal Wezero. 
The word format may:be written symbolically as 

X .XXXXXXXXXXXXX0000000000 

where X indicates wee or zero. Since the magnitude of the 
input sample is less than one, this word can always be 
treated as a 23-bit integer (excluding the sign bit) in the 
computer, and the decimal point may be inserted mentally 
whenever it is desired to determine the real value. For 
example, lt the #raction (0-5) 49 = (0.01), LSomagdeieaeco the 
fraction (0.25) 44 = (0.01).,, then the computer performs the 
addition 

010000000000000000000000 

001000000000000000000000 

011000000000000000000000 
a, = (6291456) 44- However if a decimal 


point is mentally reinserted after the sign bit, the correct 


The answer is (11x2 


answer is seen to be (0.110°**0), == (0.75) 49: Lt ‘the 24" put 


word in the answer above is input to the D/A converter, it 
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Gorrectly interprets thewvalue) as 02/5 and, the outputyis 75V. 
Thus, no conversion is necessary for data which is input via 
the A/D converter. 

The situation 1S quite different for the multiplying 
coefficients, as and b.. Lf a cock iicrventwor (0.5) 46 were 
to be entered into the program as data, the integer (base.) 


equivalent of (0-5) 49, which is zero, would be stored in 


memory. This 1S because the assembler converts all base 10 
data to base 2 integers when using integer arithmetic. The 
solution is to multiply (0-5) 49 by a. If this were done, 


then the 23-bit integer that 1s generated in the computer 
will be No times too large which is exactly what is desired. 

(0. Sales = (4194304), , = (010000000000000000000000) . 
Hence, if the base 10 integer 4194304 is entered as data the 
desired fraction is stored in memory. 

Since the as and b. for most cascade and parallel 
filter forms are within the range -2 to +2, a more convenient 
conversion is 

C, = ¢,x2°* = c, (4194304), (5-1) 
where the general notation Cs inclLudes,.aligmudstiglvangwcoef— 
ficients. The conversion in Eq. (5-1) is equivalent to di- 
viding the coefficient by two to preclude any overflow in 
the A register. After the multiplication is complete the 
answer is multiplied by two to obtain the correct numerical 


result. In Eq. (5-1) Cc. is the base 10 integer entered in 


the computer program. 


G3 


4, A Discussion of the Program 
Computer Program 2 contains the scaled coefficients 
for the second-order low-pass filter which was designed in 
Secrren Lilt Cc. The digqital™= transfer fumceron for that 
filter is 


0. 0CTA5 Sem 


— 0.067455 + 0.134910z. 
a ; 


H(z) =—aal = 
1 - 1.142980z me 4 eee Lz 


(5.525 


and the corresponding difference equation is 
CY a-Si Pe pe ty |S eo Re OS 


If the coefficients are scaled according to Eq. (5-1), 


the result is 


ay = 0.067455 Ay = 2029 27 
a, = 0.134910 A, = 565854 
a, = 0.067455 A, = Loe eT 
b, =-1.142980 By = @7294006 
b, = 0.412801 By =-1731413, 
where the signs of the B. are reversed to permit addition 


instead of subtraction.-t 


The sample values are stored in memory in the fol- 


lowing order 


DATA = x 
BATAS | =a 

n-l 
DATAt+2 = X 

n-2 
SUM mil 
SUMe 1. as Bice 9) 
SUM+2 = ve pie 


1l The reversing of the signs of the B:. in the program 
applies to all implementation schemes and will be done for 
all remaining examples without comment. 
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Then the sequence of instructions from line 22 to line 41 


in the program solves the following equation: 





sum = 2° - (DATA) +2 + 2° + (DATA+1) -2 + 3. (DATA+2) +2 
(o-4) 
+ 2. (SuM+1) +2 + 54+ (SUM+2) +2 , 


wirbeh 1s equivalent "ce Eq. (5-3). Thus;™@empucer Pregram 2 
is a Simulation of the first of the two direct filter struc- 
tures discussed in Section IV, A. 

Notice that all the coefficients were divided by two 
to permit a more general program. Also recall that, because 
integer arithmetic is being used, the most significant portion 
of a product is contained in the A register. (The entire 
product is in the A and B 48-bit register.) If the contents 
of the B register are dropped completely, a correct 24-bit 
product is still retained since the value is interpreted as 
e=meinary fraction. 

At the end of the arithmetic sequence the variables 
are shifted in memory to prepare for the next sampling inter- 
val. This is accomplished by the instructions in lines 43 
through 46. This sequence is performed five times (until the 
contents of the index register become negative), and the 
effective addresses of the load and store instructions are 
decremented by one for each iteration. The resultant shift 


in memory is shown below: 
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Sinha Order of Occurrence 


DATs last 
y 

DATA+1 Sih 
1 

DATA+2 4th 
¥ 

SUM Srd 
+ 

SUM+1 2G) 
1 

SUM+2 lst 


The current DATA+2 and SUM+2 values are lost which is of no 
consequence. The contents of SUM become correct in the next 
sample period after the STA instruction is used in line 25. 

Implementation of more complicated filters will be 
illustrated in the examples which follow. The only differ- 
ences in the program are the sequence of arithmetic opera- 
tions and the arrangement of the data and the coefficients 
in memory. 

In the beginning of this section on digital-filter 
simulation it was mentioned that it would be desirable to 
control the computer word length. This can be done conven- 
iently in the program by the use of the ETR (extract) 
instruction. Refer to line 64 of Computer Program 2. The 
ETR instruction performs a logical AND operation between the 
contents of the A register and the operand of the ETR in- 
struction. In the D/A subroutine the output value Vique'S 
loaded into the A register. Then a logical AND is performed 


between Vie and the contents of MASK which are 


(111111111111111000000000) 
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The result is that Ss uSy truncated to 15 brits waren is 

the word size used in the D/A converter. Tf an ETR instruc- 
tion is inserted after each multiplication, then all the 
arithmetic in the digital filter is truncated to the number 
of bits specified in location MASK. This means that the 
effective word length of the computer can be controlled by 
altering the contents of location MASK. Additionally the 
accuracy of the multiplying coefficients can be reduced by 
truncating the number of significant bits in the storage 
locations. This would normally be done manually since the 


Geetricients are constant in time. 
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VL... EXPERIMENTAL RESULTS OF REAL-TAME DIGRRAL—EI LTER 
SIMULATIONS 


The digital-filter simulation scheme described in the 
previous section was used successfully for the real-time 
simulation of various digital filters. The same basic pro- 
gram may be used for all filters. Different filter struc- 
tures are simulated by changing the sequence of arithmetic 
operations. Different filters of the same structure may be 
simulated by changing the filter coefficients in the program. 
For the remaining examples only the arithmetic and data- 
shifting sequences and the filter coefficients are shown in 


the referenced computer programs. 


A. A LOW-PASS FILTER 

A second-order low-pass filter with a cutoff frequency 
equal to 200 Hz and a sampling frequency equal to 2000 Hz was 
designed in Section III, C. The transfer function of that 
filter is 


dk 2 


+ 0.067455z 
2 


_ 0.067455 + 0.134910z_ 


Herz) =F ~ 
= Lea OZ + Opes 0tz 


(6=15) 


This is the same filter that was used to describe the general 
Simulation program. Computer Program 3 shows the arithmetic 
and data-shifting instructions necessary to implement this 
filter in the second-order canonical form. Notice that in 
this case additional scaling of the input data became 
necessary. This was easily accomplished using the ARSA in- 


struction in line 23. The difference equations which are 
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f. = 100 Hz 
Time = 2 msec/div 


Amplitude = 5V/div 





E. = 200 Hz 
Time = 2 msec/div 
Amplitude = 5V/div 





Figure 6-1. An Output Signal from a Low-Pass Digital Pasleer. 
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implemented by Computer Program 3 are a 


SUM _ SUM | Bl, SUM+1,, B2.,SUM+2, at 
2 ACW SUM, . Al. ,SUM+1. A2. SUM+2, . iy 
B= 5 (A) 2 + Se (6-2b) 


and the data is stored in memory as 


SUM => W 
n 
SUM+l1l = w 
n-lL 
SUM+2 = w 
n-2 
bie ~ 2a 


where WwW 1S an intermediate computational variable. 

Equations (6-2) are identical in form to the general 
second-order difference equations for the canonical form, 
Eas. (4-7). The first term on the right side of Eq. (6=2a) 
is labeled SUM only for convenience in the program. This 
location contains the current input value at the time line 
28 1s executed. Also since Y/4 1s the answer obtained from 
the: program, P036 in Fig. 5-1 must be set to 0.4 to complete 
the magnitude scaling of the output. 

‘The actual input and output waveforms for this filter are 
shown in Fig. 6-la. The ripple in the output wave 1S caused 
by the presence of the sampling frequency and its harmonics. 
These high frequencies could be eliminated by additional 
analog filtering if a smooth waveform is required. Figure 
6-lb shows the attenuation of the output at 200 Hz. The 
complete frequency response of the filter is exactly identical 


to the frequency response calculated for the transfer function. 
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Figure 6-2. The Spectrum of the Output Signal. 
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That response was piocted in Fig, 3-3. [The phase SUglaLings 
caused by the insertion of the fiiter is also evident in 
Fig. ii. 

Figure 6-2 shows the spectrum of the output waveform for 
an input signal frequency of 200 Hz. The presence of the 
spectral component at 1800 Hz (2000 Hz - 200 Hz) verifies 
that the frequency spectrum of the filtered signal is indeed 


periodic with a period equal to the sampling frequency. 


B. A BANDPASS FILTER 

A bandpass digital filter was designed in Section III, C. 
The filter characteristics included a flat passband from 
200 Hz to 500 Hz for a sampling frequency of 2000 Hz. The 
transfer function of the filter is 


A (6-3) 
4 ‘ 


DTS 1IGG20me6 oles -40,131106c" 


HD = =) 5 = 
1-1 Wow G42 © 92722 1lee 06584192 47+0 92722172 


1. Simulation of the Direct Form 
Computer Program 4 gives the abbreviated coding for 
a direct form of this filter. The difference equation, in- 


cluding the in-program scaling of line 23, 








SUM _ AO, DATA) ,, , Al, DATAr1,,, _ A2, DATAr2,,, 
ee mm oo. >. 
2s DATA+3 a DATA+4 Bl .SUM+1L 
a ws Co aa. 2 oo 2 (6-4) 
SI eu B38, sum B4., SUM+4_ . 
2 eee oe ee ce 


Figure 6-3 shows the output waveforms of this filter at 


several frequencies. The complete frequency characteristic 


We 


Be = 200 Hz 
Time = 1 msec/div 
Amplitude = 5V/div 





f. = 320 Hz 


Time = 1 msec/div 
Amplitude = 5V/div 





t. = 650 Hz 


Time = 1 msec/div 
Amplitude = 5V/div 





Figure 6-3. Several Output Signals from a Bandpass 
Digital Filter. 
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does agree with the plot of the transfer-function character- 
iSewe giver@in wae 324. 
2. Simulation of the Cascade Form 
The cascade filter form requires the factorization of 


the digital transfer function into a product of second-order 


subfilters. The result for the present example is 


=i =) 9 ( Gas) 
(0.362085-0 .722RR70 2 me Oe) 
il 


HZ) = as =) = 
(1-1.2220232 ~+0.6@3825z2 ~ )syeeoR 7 s04iz 


> * 


+0.450821z ~) 


The implementation of the cascade form involves the 
routing of an input signal through a cascade of subfilters 
(two in this case). If the same basic second-order filter 
is used for all the subfilters, then when the signal comes 
out of one subfilter it must be routed around to the input 
of the common filter as the filter coetficients are changed. 
This is the most appealing feature of the cascade form. 
High-order filters can be implemented by time multiplexing a 
low-order filter section. 

The real-time simulation scheme can be programmed to do. 
the same thing. Refer to Computer Program 5. A basic 
second-order canonical-form filter is coded with an addition- 
al indirect addressing feature. Index register 3 adds either 
One or zero (depending on which of the two subfilter iter- 


ations iS in progress) to all operands which are followed by 


the characters ",3". The data is stored in memory as follows: 
SUM = output of first subfilter = input to second subfilter 
SUM = input sample = input to first subfilter 
SUM+l = delayed sum of second subfilter 
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SUM+l1 = delayed sum of first subfilter 

SUM+2 = delayed sum of second subfilter 
SUM+2 = delayed sum of first subfilter 

Ne = final output of complete filter 
Y = ‘GelEpueNoL “first subfilter, 


and the coefficients are stored as shown in the program. The 
sequence of instructions from "REPEAT" to "BRX REPEAT" are 
‘the basic arithmetic instructions for a second-order canon- 
ical filter structure. This sequence is executed twice 
during each sample period. Indirect addressing causes the 
proper storage locations and filter coefficients to be used 
for each subfilter. In this way a time-multiplexed filter is 
Simulated. Care must be taken to insure the proper shifting 
of the data between sample intervals. Refer to lines 5l 
through 56 in the program. 

The bandpass digital filter under consideration was 
Simulated using Computer Program 5. The resultant frequency 
response was identical to both the calculated response and 
the response of the direct-form simulation. 

3. Simulation of the Parallel Form 

The parallel form of a digital filter may be chosen 
1f the paramount consideration is speed in lieu of conserva- 
tion of hardware. In this scheme the total response is the 
sum of the responses to the individual subfilters which are 
arranged in parallel. Refer to Fig. 4-4. Evidently the cal- 
culations for each subfilter can be performed simultaneously 


if no multiplexing of hardware is desired. 
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The implementation of the parallel-riiter structure 
requires the partlabefractionVexpansion of tHe digital trans-— 
fer function into a summation of second-order transfer 
functions. This operation may involve considerable labor if 
the order of the filter is high, Computer Program 6 is a 
relatively simple FORTRAN program for obtaining the partial 
fraction expansion of some rational polynomials. The multi- 
plicity of poles of the rational function cannot exceed two. 
Fortunately digital transfer functions normally do not have 
multiple poles, so the program can be used. Minor modifica- 
tions to the program would permit expansions of arbitrarily 
high-order transfer functions. 

The bandpass filter presently being considered is of 
fourth order. Accordingly, the transfer functionsgaven mm 


Eq. (6-3) must be expanded into the form 


a, (2) a, (2) 


E(z) * Bota) logs) 


H(z) = yt 


The algebra will not be given here. The results are: 


Y _ 0.131106 
a, (2) = -0.3837542 + 0.149628 
8, (z) = 2°-1.222023z + 0.603825 
a, (2) a 0.567310z2 + 0.046308 
8,(z) = 2°-0.178041z2 + 0.450821. 


Notice that the expansion was made for H(z) egual to a func- 
iO Da Oteet Zee Slee tha 7 Then, after Eq. (6-6) is written 
Lnesttondard fom (a. func .onwaik Smear the filter and program 


coefficients become 
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y = 0.131106 GAMMA = oe Ske BeRs: 
a a as. 38.34 Al, =. — 1.60053 1: 
= = 2 
Ay > 055.6 75 80 5 379471 
ao, = 0.149628 A2, = 627585 
Aan5 = 0.046308 a 194230 
= gs Bl = 512558 
11 Lu 222028 1 6 
Di5 = -0.178041 aa 746758 
boy = 0.603825 B2, = -2532Z020 
Do5 = 0.450821 a —eSiOroe'O 


Computer Program 7 is the simulation of the parallel form 
of the bandpass digital filter. The program is identical to 
the program for the cascade form except that the outputs of 
each individual subfilter are stored separately and summed 
together at the end. The simultaneous operation of the 
subfilters obviously cannot be simulated because the com- 
puter's single arithmetic unit must be time-multiplexed. The 
correct results for the parallel filter are obtained although 
the increase in operating speed is not. In fact, the parallel 
filter simulation takes slightly longer than the cascade 
simulation because of this additional summation sequence 
(lines 49 through 53) which is required at the end. 

The results of the real-time simulation of the bandpass 
filter using Computer Program 7 were essentially identical 


to those of the previous two simulations. 


4. Experimental Observations of Quantization Effects 


The multiplying coefficients were calculated to six- 
decimal-place accuracy in the design of the bandpass digital 
filter. This accuracy 1s retained by the 23-bit magnitude 


representation in the XDS 9300 computer. If the least 
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significant bits of the coefficient data words are system- 
atically truncatedy some Giisight wnto the Leeeets of 
coefficient quantization can be obtained. A software package 
is available with the XDS»9300 computer which» permits alter- 
ing the contents of memory locations vila manual input to the 
teletypewriter. 

The coefficients for each of these bandpass-filter 
configurations were truncated in steps of single octal 
integers with the proper rounding being performed. Generally 
the effect of the less accurate coefficients was to gradually 
alter the frequency-response characteristic until it no 
longer satisfied the filter specifications. This occurred 
when the program coefficients were reduced to three signif- 
icant octal digits 1n each case. Thus, eight magnitude bits 
or third-decimal-place accuracy was required for proper 
operation of the bandpass digital filter. 

The results of this example do not permit any general 
conclusions to be made regarding the choice of one filter 
form over another. It is likely that the arguments against 
the direct form from the point of view of coefficient accuracy 
which retinas in the literature [{l, 4, and 5] have greater 
Significance in higher-order filters. Also the problem of 
filter stability was not encountered in this example, since 
the poles of the digital transfer function were located well 
Within the unit circle in the z-plahe. 

MOGrrreatron OL the a oe length 1s accom- 


plished in the simulation program by inserting the ETR 
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instruction after each multiplication and altering the 
contents of location MASK as desired. Although the A/D 
converter provides only 13 bits of accuracy in the input 
data, all calculations are performed with 24-bit accuracy 
unless altered by the ETR instruction. Accordingly the 
effects of round-off error in a digital filter can be ob- 
served as the computer word length is reduced. 

The contents of location MASK were decremented from 
15 bits for each filter form. The point at which quantiza- 
tion noise reduces the response of the digital filter to an 
unsatisfactory state must be judged in light of the filter 
application. For this example eight bits and seven bits 
were required for the direct form and the parallel form 
respectively to retain basic signal definition on an oscil- 
loscope. The response of the cascade filter was unsatisfac- 
tory even when the computer word was truncated to 15 bits. 
Figure 6-4 shows the output waveform for a 300 Hz signal 
using 15-bit arithmetic in the cascade and parallel filters. 
Figure 6-5 shows the relative quantization noise in the 
spectra of the two signals. 

It appears that the cascade Eee structure, with 
its serial form of signal processing, requires considerably 
greater accuracy in arithmetic computations than either of 
the other two structures. Hence, the significant advantage 
of using time-multiplexed circuitry is somewhat offset by 


the increased word length requirement. 
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f. = 300 Hz 


Cascade Form 


Figure 6-4. 


f. = 300 Hz 


Parallel Form 





Quantization Effects in Digital Filters 
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C. A NOISE FILTER 

As a final example consider the problem of generating a 
noise signal whose power spectrum is of uniform height within 
a specified frequency range, but 1s significantly attenuated 
outside this frequency range. Such a noise signal has 
application in various signal-analysis and measurement schemes 
for communications systems, radar systems, etc. 

If a continuous-wave carrier signal is frequency-modulated 
by a gausSian noise source, the resultant signal spectrum is 
Belssian Uimemtaeemancd Cemeercason the carrier frequency. The 
desired noise spectrum could be obtained if this gaussian- 
shaped spectrum could be filtered by a bandpass filter whose 
magni tude-ver sus-frequency characteristic varied as the 
reciprocal of a gaussian curve within the passband. 

Such a problem in spectral shaping can be solved using a 
digital filter. The bandpass digital filter which was pres- 
ented in Section III, C, as an example of matched z-transform 
design is offered as an approximate solution to the problem. 

The gaussian spectrum was generated by applying a low- 
frequency noise source (the amplitude distribution of the 
noise voltage is gaussian) to a voltage-controlled oscil- 
lator. The carrier frequency was chosen to be 450 Hz and 
the rms voltage of the noise was 1.0 volts. Applying a con- 
stant 1-0 volts to the VCO displaced the carrier frequency 
200 Hz. Therefore the gaussian spectrum centered at 450 Hz 


has a standard deviation of 200 Haz. 
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f. = 300 Hz 


Parallel Form 





ff. = 300 Hz 


Cascade Form 





Figure 6-5. Quantization Effects in Signal Spectrums. 
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The bandpass second-order Chebyshev filter was 
designed for a 12-db gain between the peaks and the valley. 
This corresponded to a displacement of 1.6 standard deviations 
from the mean on the normal curve; or, in terms of frequency, 
for the given hardware it corresponds to the peak frequencies 
being separated from the carrier by 320 Hz. Thus, the ideal 
frequency characteristic for the specified conditions is 


shown in Fig. 6-6. The actual solution obtained by standard 


10 db 


@ 


30H Zeee25 Onli z, 450 Hz 6505He “770 Hz 


Figure 6-6. An Ideal Noise Filter. 


analog-filter design methods and a subsequent matched z- 
transformation was plotted in Fig. 3-7. Notice that the s- 
domain solution in Fig. 3-6, which is itself distorted by the 
low-pass-to-bandpass transformation, suffers further distortion 
by the matched z-transformation. The final result is, however, 
a fair approximation of the desired frequency characteristic. 
The gaussian spectrum, centered at 450 Hz, is shown in 
Fig. 6-7a and the filtered spectrum is shown in Fig. 6-7b. 
Notice that the gain factor in the digital filter was adjusted 


for 6db gain at the peak frequencies and a 6-db attenuation 


103 





a. Gaussian Input Spectrum 





b. Filtered Output Spectrum 


Figure 6-7. Output Signal and Spectrum from a Digital 
Noise Filter, 
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at the center frequency. Murther refimement inthis £1i¢er 
design could be accomplished using optimization techniques 


in the s domain. 


VII. SUMMARY AND CONCLUSIONS 


The subject of digital signal processing continues to 
receive ever increaSing attention. The basic appeal of 
digital circuitry is still lower cost, higher reliability 
and faster operation. Specifically, real-time digital filters 
now appear to be on the threshold of becoming a commercial 
reality. To a large extent the impetus for the theoretical 
development of the digital filter has been the promise of 
special-purpose computational hardware using LSI technology. 
It appears extremely likely that digital filters will find 
significant real-time applications in the not too distant 
future. 

An attempt has been made in this thesis to assemble and 
organize much of the current theory of digital filters, and 
to present in detail those aspects relating to the real-time 
implementation of digital filters. Furthermore, the proposed 
hybrid-computer simulation scheme provides a vehicle to 
process signals in real time and to observe the theory of 


digital filters in action. 


A. COMMENTS ON DESIGN METHOD AND IMPLEMENTATION TECHNIQUES 
Most of the theory of the z-transform and its application 


to discrete-time signals is not new. However, the development 
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of the digital filter has brought about new variations of 
that transform. For example, discussions of the bilinear 

and the matched z-transforms are not found in the older texts 
on sampled-data systems. Also the very popular digital- 
Filter design technique of using s-plane continuous-filter 
synthesis has generated new and extremely useful "special 
purpose" transformations. One of these transformations, the 
bandpass version cf the bilinear z-transform, obviates the 
conventional low-pass-to-bandpass frequency transformations 
in the s-plane. 

The usefulness of all the various transformations lies in 
their ability to be easily implemented and to produce accu- 
rate results. The use of large computer programs to synthe- 
Size digital transfer functions from continuous transfer 
functions adds versatility to the indirect design methods. 

The three general implementation schemes which have been 
discussed in this paper all have useful applications. The 
direct form (including both the straight and canonical 
versions) is the simplest to use and gives very satisfactory 
results for relatively low-order filters. The cascade form 
represents a high-order filter as a cascade of low-order 
subfilters. This scheme is particularly adapted to realiza- 
tion in special-purpose hardware because of the "basic filter 
section" concept. A second-order filter section is time- 
multiplexed to permit the cascading of a single signal 
through the filter section n times. In this way a 2nth-order 


filter can be realized. 
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Time multiplexing is the key to the economical use of 
hardware. For example, the serial multiplier which is 
presented in Appendix A would seem to be a disproportionately 
complex circuit were it not for the built-in ability to 
handle multiplexed input signals and coefficients. 

The final filter structwre is the parallel form. Thais 
form could provide the fastest signal processing if separate 
hardware were available for each subfilter. On the other 
hand, the parallel structure can also make use of time- 
multiplexed hardware. This implementation scheme retains 
the compact and economical construction of the cascade form, 
but could operate with a shorter word length. This is be- 
cause the cascade form inherently requires greater accuracy 


in the arithmetic unit. 


B. AN APPRAISAL OF THE REAL-TIME DIGITAL-FILTER SIMULATION 
The hybrid-computer program which was presented in this 
paper represents a first step toward the real-time processing 
of signals by digital filters. Theoretically any digital 
filter that is designed for real-time application could be 
simulated by this general scheme. (Real-time digital filters 
are considered to be of the recursive type in this context.) 
In the author's view there are two basic limitations to 
the digital-filter simulation scheme. The first is the 24- 
bit word length of the XDS 9300 computer. The coefficient 
and arithmetic accuracy required in high-order filters may 
not be attainable using a =24-bit computer word. A seluceen 
to this problem would be to rewrite the simulation programs 


for double-precision arithmetic. 
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The second limitation is that there 1S an upper iim. c 
on the allowable sampling frequency imposed by the cycle 
time of the XDS 9300 computer. Execution of a second-order 
difference equation (one basic filter section) in Computer 
Program 7 requires approximately 200 microseconds of computer 
time. This figure varies slightly depending on the arithme- 
tic sequence being used. This minimum execution time means 
that a fourth-order filter cannot be simulated using a 
sample frequency higher than 2500 Hz, or that an eighth- 
order filter cannot be simulated uSing a sample frequency 
hagmemethan 1250 Hz, etc. From these figures 2t Ts "evideme 
that experimentation must be limited to the sub-audio range 
of frequencies. Indeed, the same general problem applies 
to all digital-filter systems. 

A partial solution to the problem at the Naval Post- 
graduate School is to write a more efficient program for 
the XDS 9300 computer. However, the future for the general 
Simulation scheme and for digital filters as a whole is 
much brighter. The next generation of digital computers will 
make extensive use of LSI technology. Computers using mon- 
OlmeitewcHreurts will have addition times in phemoner of 
50 nanoseconds vice 1.75 microseconds for the XDS 9300 com- 
puter. This reduction of almost two orders of magnitude in 
computing speed will mean an expansion of the digital- 
Filtering spectrum to include the complete audio range and 
beyond. There is every reason to believe that special- 


purpose hardware will benefit from similar advances. 
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APPENDIX A 


A BENARY SERIAL MULTIPETERWBOR, APPLICATIONS IN DIGITAL FLLTERS 


alli INTRODUCTION 

Real-time implementation of digital filters involves the 
construction of a special-purpose digital computer. The 
Brithmetic unitwotethbacecomputer, whicheperformss tyeamuiter.= 
plication and addition of digital signals, presents some 
interesting challenges in logical design. One part of the 
arithmetic unit, the multiplier, is the subject of this 
Appendix. 

Specifically the computational aspects of an arithmetic 
unit which are peculiar to digital filters are discussed. 
From these considerations some general specifications for a 
multiplier are determined. Finally one solution to the prob- 


lem 1S presented using a binary serial multiplier. 


2. COMPUTATIONAL CONSIDERATIONS IN DIGITAL -FILTERS 

Most digital-filter transfer functions permit simultane- 
ous arithmetic operations during a given Nyquist interval. 
If the number of operations to be performed is large, economy 
can be achieved using serial arithmetic and by sharing the 
adders and multipliers. Time sharing in the filter can be 
accomplished with multiplexing circuitry. In addition to a 
Significant simplication of the hardware, serial arithmetic 
provides increased modularity and flexibility WH thie. difeital 


circuit configurations. Also the processing rate is limited 
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only by circuit speed and not by carry-propagation times in 
the adders and multipliers. Finally with serial arithmetic 
sample delays are realized simply as single-input single- 
output shift registers. 

The two's complement representation of binary numbers is 
also appropriate because addition can proceed starting with 
the least significant bit with no advance knowledge of the 
Signs and with no later correction to sums being necessary. 


For example, the algebraic summation of the two signed num- 


bers, 1.101 and 0.010, can be accomplished as follows: 


conversion to 
two's complement 


eee Ou O 10a, 
Subtract 0.010 LS | SGnL@ Add 
0 . Ome: lion 


Addition using the two's complement form gives the correct 
Signed result with no concern about signs or carries. tin 
the multiplier, however, it 1s simpler to treat the sign 
separately and multiply only positive numbers. After the 
positive product is obtained the correct sign is reinserted 
iD Giemppopere Olt Location. 

Multiplication in digital filters is subject to a very 
severe restriction. Normally the multiplication of an N-bit 
number by a K-bit number results in a product (N+K) bits 
long. This situation cannot be tolerated ina serial system 
because the N-bit words (representing successive samples) are 
adjacent to each other in time. Consequently there is no 
room for the generation of a product in excess of N bits. 


If an N-bit sample is multiplied by a K-bit coefficient, the 


15516, 


product must be truncated and rounded to N bits. The solu- 
tion to this problem is to multiply one bit at a time using 
partial sums. It can be shown that the truncation of each 
partial sum and rounding the final partial sum will produce 
the “same result as trmiiicatang’ and» round iwgeasiul) fine 
proauct 


EXAMPLE 


Let x the sample value = 0.00101, where the last bit 


in time is the sign bit (0 always). 
a = the multiplying coefficrvent’="1. 0190 watthwne 
Sileneme came 
In this example N=6 and K=5. Then multiplication in the 


normal fashion can be compared to partial sum multiplication 


as follows: 


NORMAL PARTIAL SUM 
at), 0c aoe Ge, Ome 
1 err io 171 
Ona) OmOmom0 000008 
0 Ome 0 1 0.0 0. te 
0 Io “Om em 7 0% OM O'R 
Ou OocOui sunt ey ee 
go UT Um or osur rt ar 
O70 OTL TOME gD OO de 0.08 
1 — 
0-00 IT woelo nr 
000K pale 
ik 
Sa Om laaL 


The results are identically correct. 

If we require the magnitude of the sample to be less than 
one (1.¢e., lsseahe< Jae Oi), and the magnitude of the multiplying 
coefficient to be less than two (i.e., |a|<2.0), then the 
fixed-point multiplication of the positive coefficient times 
the positive sample will always be of the form shown in the 


example. That is, the product is always positive and the 


ae i. 


decimaimpoint follows sthemsignebia teguweteaqsms vedid in, thie 
input «er the multiplrers eit is» important tomoterthat the 
stom bit of tthescoefficient»isenot ancbuded in ithe mudéa- 
plication. The restrictions on magnitude size are not 
unrealistic because the sample can easily be scaled. Also 
the multiplying coefficients for typical cascade and parallel 
filter forms are in the range -2.0 to 2.0. Any overflow in 
the partial sums will cause an error. This problem can be 
eliminated by further scaling and will not be discussed. 

In view of the above discussion the general requirements 
for the multiplier are the following: 


l. Accept an uncomplemented sample value from a serial 
delay device. This sample will be in binary form, N bits 
long, with the last bit in time representing the sign. The 
remainder of the word will be a fractional value less than 


one. 


2. Accépt! each single bit of a K-bit multiplying coef- 
ficient separately for use as a constant bit multiplier. The 
Sign bit is used only to determine the sign of the product 


and the magnitude of the coefficient is less than 2.0. 


3. Remove the sign bit, from the sample value and 


x 
NN” 
multiply it by the coefficient sign bit, ae modulo 2 


(exclusive-OR) to determine the sign of the product. 


4. Multiply the remaining positive sample by the positive 
multiplying coefficient using a serial multiplication scheme. 
The partial products must be truncated to N bits during éach 
step and the final N-bit product is rounded up by adding in 
the last truncated bit. The sign bit of this product must 


always be a zero (positive). , 


9. Reinsert the proper sign bit into the final product. The 
decimal point will be positioned immediately after the sign bit 


in the product. 


lite2 


6. Convert the signed product to the two's complement form 


for insertion into the following adder. 


3. THE MULTIPLOSR 

A binary serial multiplier which will satisfy the above 
requirements is shown in block diagram form in Fig. A-l. An 
explanation of the multiplier follows: (Refer to correspond- 
ing step numbers on the diagram). 


1. The notation of Fig. A-l is based on a sample value N 
bits long inéluding sian aud a multiplying coeffictent K bases 
long not incl¥dingmeagn. | "Tame ee is defined as the time the 
first bit of the sample, X,, becomes available at the input. 
This is the least significant bit of the sample word. Then 
succeeding times become t,, to,ctety_y (which is the time the 


last bit of the sample, becomes available at the same 


Xe 
Bount) « 

2. First the sign bit of the sample is removed and that 
bit 1s converted to zero (positive) in the sample word. This 
operation normally requires an N-bit delay but the delay will 
be disregarded for sake of clarity in the discussion which 


follows. 


3. The removed sign bit, is multiplied modulo 2 with 


My 


the sign bit of the coefficient, a The resulting product 


Sign is stored or delayed for Kona ani 
4. At time to the least significant bit of the sample, Xa 
is multiplied by the least significant bit of the coefficient, 
ay: This product tsmene first bit of *He first” partial semm:- 
The inputs to the multiplier bit sections labeled vr (7=0Rl , 
P*-K-2) are timing’ sigmals witch accomplish) the trunedeton. 
At time to. ©, would be low (the other vs being high) so the 

output of partial "sam IVMae time to would be gated to zero. 
At time t, r, goes low so the output of partial sum 2 at ty 


is zero. In this fashion the signal at the input to the 


JAS 


(laa ements BatwiSectaions 





Multiplier Bj a ion 


miltep 


DATA 





Old Partial Su 
New pearl a Sum 


a. ‘ 
aL. v5 


Note: = is low at bit time j. The line is high all other 


times. 
Figure A-l. A Binary Serial Multiplier. 
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multiplier follows the sequence xX, (ty), Xo (ty) ,oe ey (ty_)) 
while the signal at the output follows the sequence O(to), 


illustrate the scheme more clearly consider the following 


example: 
EXAMP LE 
Let x=0.010 N=4 
a=1.011 K=4 
Then x = Coa 0 
a= Ome) 
© O01 0 Ist "partial stm 
Ome 0 
Oo Mil Z 2nd partial sum 
00 0 0 
00027 3rd partial sum 
020530 
1 Last truncared ome 
y = 0-041 4th partial sum 
Cg C5 eum ec oe cy 


(K-1) bit delay in product. 


It 1S easily seen that the preceding sample, occurring at 
the input times t_yar tla, t t_* would have: predteed a 


-2' ~“=1 
correct product at the output at times tlie as tis t.- 


5. In the last multiplier bit section, the bit which is 


truncated at bit “Gime tei is»ted backethrough wie rama 


loop to accomplish the rounded product. 


6. Finally the sign bit (zero in all cases except over- 
flow) is changed to the correct value and the correctly 
Signed product is then converted to two's complement form for 


insertion into an adder. 

It is seen from Fig. A-1l and the above discussion that 
this serial multiplier has an inherent delay of (K-1) bits, 
not counting any delay encountered in the sign bit removal. 
This additional delay must be deducted from the prescribed 
delay in the multiplication branch for proper operation of 
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